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ABOUT  THE  STUDY  IN  THE  TIME  DOMAIN 
OF  JUNCTIONS  BETWEEN  THIN  WIRES 


Amelia  Rubio  Bretones 
Alfonso  Salinas  Extremera 
Rafael  Gomez  Martin 
Jesus  Fornieles  Callejon 

Dpto.  Flsica  Aplicada 
Facultad  de  Ciencias 
Universidad  de  Granada 
18071 -Granada  (SPAIN) 


ABSTRACT 

DOTIGl  is  a  computer  code  developed  for  the  study  of  the  interaction 
of  arbitrary  electromagnetic  signals  with  thin-wire  structures,  in  the  time 
domain.  It  calculates  the  current  distribution  induced  on  the  structure  by 
solving  the  electric  field  integral  equation  using  the  moment  method.  The 
numerical  procedure  used  to  develop  the  program  and  different  possibilities  for 
treating  junctions  are  briefly  described. 

To  obtain  an  accurate  solution  for  the  current  induced  on  the  thin-wire 
structures  it  is  very  important  to  pay  attention  to  the  zones  at  which  the  wires 
intersect.  Thus,  different  junction  treatments  were  tested  for  several  simple 
structures.  Following  some  convergence  criteria  the  current  distributions  were 
compared  to  a  reference  solution  and  also,  by  way  of  Fourier  transform,  with 
results  obtained  using  some  well  known  frequency-domain  codes. 


1.  INTRODUCTION 


Recent  developments  in  high-resolution  radar  and  electromagnetic-pulse  technology 
(EMP)  have  generated  interest  in  their  interaction  with  structures.  This  can  be  studied  either 
by  using  time-harmonic  analysis  followed  by  Fourier  inversion  or  solving  directly  in  the  time 
domain  using  a  time-step  algorithm  [1].  The  time  domain  approach  not  only  offers  computational 
advantages  but  sheds  light  on  the  problem  in  a  way  that  is  not  possible  using  frequency-domain 
techniques  [2]. 

The  first  step  for  such  an  analysis  is  to  calculate  the  electric  current  induced  on  the 
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surface  of  the  structures.  As  they  are  normally  very  complex,  a  common  technique  is  to 
approximate  the  surface  geometry  by  a  set  of  joined  wires.  To  obtain  an  accurate  numerical 
solution  for  the  structure  current  distribution,  particular  attention  must  be  paid  to  the  zones  at 
which  the  wires  are  joined  together  [3]. 

In  this  paper  different  treatments  for  the  study  of  the  junctions  are  presented  in  order 
to  compare  how  stable  and  accurate  the  solution  is  obtained.  The  treatments  were  applied  to: 
a  cross  of  wires,  a  stepped-radius  wire,  a  net  (Fig.  13)  and  a  dodecahedron  (Fig.  16),  when 
illuminated  by  a  transient  electromagnetic  field.  The  currents  induced  on  the  wires  were 
calculated  in  the  time  domain  by  solving  the  electric  field  integral  equation  (EFIE)  using  the 
moment  method.  The  results  obtained  using  the  different  junction  treatments  were  compared, 
by  Fourier  transform,  with  those  obtained  from  the  computer  code  NEC-2  [4],  working  in  the 
frequency  domain,  and  for  the  dodecahedron  with  the  AWAS  code  [5]  as  well. 

2.  NUMERICAL  METHOD 


DOTIGl  computes,  in  the  time  domain,  the  currents  induced  on  a  structure  modelled  by 
interconnected  wires  when  it  is  excited  by  an  arbitrary  electromagnetic  signal.  The  method 
employed  consists  of  resolving  the  electric  field  integral  equation  (EFIE)  using  the  moment 
method.  This  equation  for  a  thin  wire  is  {6],[7]: 


S-eUs,  t)  =—^-( 

4iieo  Jc(s>)  c^R  ot  cR^  os 


(1) 


R 

R^ 


g{s',  tO  ]  ds' 


where  s  and  s’  are  tangent  vectors  to  the  wire  axis  of  contour  C(s’)  at  position  s(r)=s  and 
s(r’)=s’.  I(s’,t’)  and  q(s’,t’)  are  the  unknown  current  and  charge  distributions  at  source  point  s’ 
at  retarded  time  t’=  t-R/c  and  E‘  is  the  field  applied  to  the  observation  point  R  =  |  r-F  |  (See 
Fig.  1).  The  charge  q(s’,t’)  can  be  expressed  in  terms  of  I(s’,t’)  by  the  equation  of  continuity. 

Using  the  point  matching  form  of  the  moment  method,  integral  eq.  (1)  was  transformed 
into  the  group  of  equations  (in  matrix  notation)  [8], [9] 


(2) 


which  allows  us  to  calculate  the  current  Ij  at  time  tj  from  the  elements  E‘j  of  the  tangential 
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electric  field  scattered  by  currents  of  previous  times  and  the  elements  E'j  of  the  tangential  applied 
electric  field  at  the  observation  point  and  at  time  step  tj.  Matrix  Z  is  a  matrix  of  interaction,  the 
elements  of  which  are  time-independent.  They  depend  on  the  geometry  and  on  the 
electromagnetic  characteristics  of  the  structure  alone  [10]. 


The  base  function  used  was  a  nine-point  lagrangian  interpolation,  which  in  the  case  of 
a  single  wire  can  be  written 


t') 


+1  +1 

■E  E 


{l,Ta) 


(3) 


where  Ii,j(s’,t’)  is  the  value  of  the  intensity  current  at  position  s’  and  time  t’  inside  one  of  the 
space-temporal  intervals  i-j  into  which  the  wire  was  divided.  are  the  coefficients  of  the 

interpolation  that  depend  on  s’  and  t’,  and  Ij+y+m  are  the  values  of  the  current  at  the  space- 
temporal  interval  i-t-l-j-f-m  nearest  to  the  i-j  interval  [8],[10]-[11]. 

Fig.  2  shows  how  the  space  interpolation  is  carried  out  for  an  isolated  wire.  The  current 
at  a  point  s’  in  the  i-segment  is  interpolated  to  the  values  of  the  currents  at  the  centers  of  its 
neighbouring  spatial  intervals  (i-1,  i  and  i+1)  (fig-  2a).  The  current  on  the  end  segments  is 
interpolated  to  zero  just  at  the  wire  ends  (fig.  2b)  [12]. 
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Fig.2  Spatial  interpolation  for  a  single  wire. 

a)  Not  end  segments,  b)  End  segments. 


3.  DESCRIPTION  OF  THE  DIFFERENT  MULTIPLE-JUNCTION  TREATMENTS. 

When  dealing  with  structures  that  have  several  interconnected  wires,  the  assumption  that 
the  values  of  the  currents  at  the  ends  meeting  at  the  junction  are  zero  is  invalid,  so  the  current 
on  the  segments  closest  to  the  junction  point  cannot  be  interpolated  to  zero  and  the  numerical 
method  described  in  section  2  cannot  be  applied  [12]. 

For  treating  junctions  we  have  modified  the  method,  trying  several  schemes  which  are 
summarized  in  Fig.  3;  the  main  differences  between  them  are  the  way  of  discretizing  the  wires 
and  how  the  interpolation  function  is  handled  near  the  junction  [13]: 

Tl:  Overlapping-segments  treatment: 

The  first  method  employed  was  the  translation  to  the  time  domain  of  a  technique  that 
had  been  used  successfully  in  the  frequency  domain  [12].  The  method  consists  of  replacing  the 
actual  configuration  of  N  wires  that  meet  at  the  junction  (see  Fig.  3a)  by  a  set  of  electrically 
unjoined  but  overlapping  wires  (Fig.  3b).  It  is  easy  to  prove  that  the  way  chosen  to  represent  the 
overlapping  wires  satisfies  ICirchoff  s  first  law  at  the  junction.  This  way  is  that  when  N  wires  meet 
at  any  junction,  wire  i  overlaps  one  segment  onto  wire  i-1,  except  that  wire  1  is  not  overlapped 
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onto  wire  N.  A  wire  is  then  considered  to  end  not  at  a  junction  but  at  the  end  of  the  overlap, 
where  the  current  is  set  to  zero.  This  way  each  wire  can  be  treated  as  if  it  were  isolated. 

The  interpolation  of  the  current  on  a  segment  will  be  to  the  sum  of  the  current  on  that 
segment  and  the  extra  overlapping  segment  if  it  exists,  as  shown  by  the  arrows  in  the  diagram 
(Fig.  3b). 

As  matrix  Z  in  eq.  (2)  depends  only  upon  the  geometry  of  the  structure,  the  application 
of  the  moment  method  to  a  segment  and  its  overlapping  segment  leads  to  exactly  the  same 
equation.  For  a  junction  of  N  wires  the  system  of  equations  (2)  has  N-1  identical  equations  (the 
number  of  overlaps).  These  extra  equations  were  removed  and  substituted  by  the  Wu  and  King 
condition  for  the  charge  at  the  junction  [14],  which  for  thin  wires  can  be  expressed  as  the  match 
of  the  derivative  of  the  current  over  space  at  each  wire  at  the  junction.  Therefore,  a  linear  system 
is  obtained  with  the  same  number  of  unknowns  as  equations. 

T2.  Half-segments  treatment  (Fig.  3c): 

Each  wire  was  discretized,  locating  at  each  end,  segments  of  half  the  length  of  the  other 
segments  of  the  structure.  The  current  on  those  half  segments  is  interpolated  backwards  to  the 
currents  on  segments  of  the  same  wire,  as  shown  in  Figure  3c. 

The  Kirchoff  s  first  law  and  the  Wu  and  King  condition  described  in  the  T1  treatment  are 
imposed  adding  extra  equations  to  the  ones  that  come  from  the  moment  method.  In  this  case  the 
method  leads  to  a  linear  system  of  equations  with  one  equation  more  than  the  number  of 
unknowns.  The  system  is  solved  by  multiplying  eq.(2)  by  Z  transpose,  which  minimizes  the  square 
error. 

T3.  Sum  treatment  (Fig.  3d): 

The  multiple  junction  of  N  wires  is  treated  by  interpolating  the  current  from  one  segment 
across  the  junction  to  the  center  of  a  "ghost  segment"  with  current  equal  to  the  sum  of  the  N-1 
currents  of  the  other  segments  at  the  junction. 
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It  can  be  seen  that  the  Wu  and  King  condition  is  implicit  in  this  junction  treatment,  but 
not  the  Kirchoff  s  first  law,  which  has  been  imposed  by  adding  an  extra  equation  to  eq.  (2)  and 
solving  as  described  in  T2. 

4.  RESULTS 


The  different  junction  treatments  described  in  section  3  were  applied  to  several  simple 
structures  in  order  to  study  the  accuracy  of  the  solution  obtained  and  the  efficiency  of  the 
method  employed. 


The  first  structure  analyzed  was  a  cross-wire  structure  (see  Fig.  4).  The  currents  induced 
on  it  were  calculated  when  it  was  illuminated  by  a  gaussian  pulse  of  the  form: 


.B=exp  ( -g^  ( t- 2 )  f 


(4) 


with  g=5*10’ s'^  and  t„„=4.28*10“’ s.  The  radius  of  all  the  wires  was  0.00222m. 

In  Fig.  5  the  current  distribution  is  presented  versus  time  at  a  point  on  the  first  wire 
situated  4.6  mm  from  the  junction  point  (Point  O). 


In  order  to  compare  how  convergent  the  different  junction  treatments  are  with  respect 
to  the  number  of  segments  in  which  the  structure  is  divided,  two  error  criteria  were  calculated: 
the  root  mean  square  current  Ifj,s  and  the  normalized  root  mean  square  current  error  e^s  [15]- 
These  variables  are  defined  by: 


^ims  O 


s,  Ht 


N^Nf 


t  i=i 


w/  Ki 


where  is  the  number  of  segments  into  which  the  structure  was  divided  and  Nt  the  number  of 
time  intervals  calculated  in  each  case.  N/,  N/  and  Pg  correspond  to  a  reference  solution  that  was 
obtained  dividing  the  wires  into  a  sufficiently  large  number  of  segments. 
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K 


wire  3 

Fig.4  Cross-wire  structure.  p=0.11  m,  a  =0.00444  m. 

Length  of  the  wires:  (wire  1=  wire  2=  wire  4=  p; 
wire3=  2p). 


Figs.  6  and  7  show  respectively  I^ns  and  e^ns  versus  the  number  of  segments.  The  rapidity 
with  which  a  technique  leads  to  a  stable  solution  can  be  estimated  from  the  plot.  It  is 
desirable  that  the  curves  should  approach  a  constant  as  soon  as  possible.  As  can  be  seen  in  Fig. 
6,  the  Inns  obtained  using  T1  and  T3  approach  a  constant  more  quickly  than  T2.  The  e^^^s  plot  in 
Fig.  7  also  shows  that  T1  and  T3  behave  better  than  T2,  that  is,  these  approaches  need  fewer 
segments  than  T2  to  arrive  to  solution  with  ^  0  . 

0.4  - 

0.3  ‘ 

0.2- 

I 

i  0.1  - 

I- 
Z 
UJ 

on  0  - 

z> 
o 

-0.1  “ 

-0.2- 

-0.3- 

0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6 

TIME  (s) 

-  OVERLAPPED  -  HALF  SEGMENT  .  SUM 

Fig.  5.  Current  distribution  versus  time  at  point  O  on  the  cross  wire  structure. 

Overlapped  and  Sum  give  the  same  results. 
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60 

o  SUM 


Fig.  6  Root  mean  square  current  versus  the  number  of  segments  into  which  the 
cross-wire  structure  was  divided. 


Fig.7  Normalized  root  mean  square  current  error  versus  the  number  of 
segments  into  which  the  cross  structure  was  divided. 

In  order  to  compare  our  time-domain  results  with  the  frequency  domain  results  given  by 
the  computer  code  NEC-2  [4]  we  Fourier  transformed  them  and  divided  by  the  transform  of  the 
incident  field.  Fig.  8  shows  the  Fourier  transformation  of  the  current  at  a  point  4.6  mm  from  the 
junction  versus  frequency.  The  results  obtained  with  T1  and  T3  agree  closely  with  those  obtained 
from  NEC. 

The  same  analysis  was  made  for  the  case  of  the  stepped-radius  wire  (considering  it  as  a 
junction  of  two  wires  with  different  radii)  drawn  in  Fig.  9,  when  illuminated  by  a  gaussian  pulse 
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CURRENT  (mA).  IMAGINARY  PART 


Fig.8 


FRECUENCy  (GHz.) 

OVERLAPPED  — -  HALF-SEGMENT  .  SUM  - NEC 


FRECUENCy  (GHz.) 

OVERLAPPED  -  HALF-SEGMENT  .  SUM  -  NEC 


Current  at  point  O  on  the  wire-cross  structure  versus  frequency,  a)  Real  part,  b) 
Imaginary  part.  (Overlapped  and  Sum  give  the  same  results). 


with  g=5*10’  s  >  and  t„aj=9.2*10-i“  s.  Figs.  10, 11  and  12  show  respectively  the  current  versus  time 
at  point  O,  the  e,ms  error  versus  the  number  of  segments  and  the  comparison  of  the  Fourier 


transform  of  the  time-domain  results  with  the  ones  obtained  by  NEC.  As  can  be  seen,  in  this  case 
all  the  treatments  behaved  similarly. 


T 


Fig.  9.  Stepped-radius  wire  geometry.  p=  0.3m,  q=  0.2m, 
a=  0.002m  and  b=  0.004m. 


TIME  (s)  *10E-8 

-  OVERLAPPED  -  HALF  SEGMENT  .  SUM 


Fig.lO  Current  distribution  versus  time  at  point  O  on  the  stepped-radius  wire.  (Overlapped 
and  Sum  give  the  same  results). 


The  following  structures  studied  were: 

1. -  The  net  of  twelve  wires  represented  in  Fig.  13. 

2. -  A  dodecahedron  structure  (Fig.  16). 

In  both  cases  the  electric  field  incident  was  proportional  to  the  temporal  derivative  of  eq.  (4) 
polarized  in  the  x-direction  and  propagating  in  the  z-direction: 
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E^=2g(  t- exp  ( -g^  ( 2)  j? 


(7) 


with  g=l*10’ s*  and  tm„=4.6*10’ s  for  the  net  and  g=8*10*  s'*  and  t„„=5.7*10'®  s  for  the 
dodecahedron.  This  incident  field  was  chosen  because  it  has  not  zero-frequency  component  and, 
therefore,  does  not  excite  the  circulating  current  that  would  exist  on  structures  with  closed-circuit 
parts  [10]. 


NUMBER  OF  SEGMENTS 

•  OVERLAPPED  +  HALF  SEGMENT  SUM 


Fig.  1 1  Normalized  root  mean  square  current  error  versus  the  number  of  segments  into  which 
the  stepped-radius  wire  was  divided.  (Overlapped  and  Sum  give  the  same  results) 

As  the  T2  treatment  was  found  to  give  less  accurate  results  when  compared  with  NEC 
and  to  need  more  segments  to  arrive  at  a  stabilized  solution,  for  these  two  more  complex 
structures  a  comparison  was  only  made  between  the  T1  and  T3  treatments. 

In  Figs.  14  and  17  the  current  distribution  is  represented  versus  time  at  point  O  on  both 
structures,  and  the  magnitude  of  the  Fourier  transformation  of  these  figures  was  compared  with 
NEC  and  plotted  in  Figs.  15  and  18.  The  results  for  the  dodecahedron  structure  were  also 
compared  with  the  ones  obtained  from  the  computer  code  AW  AS  [5].  In  both  cases  the  results 
arrived  at  with  the  T1  treatment  (overlapping  segments)  were  found  to  agree  with  NEC  and 
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-  OVERLAPPED  .  SUM 


Fig.  14  Current  distribution  versus  time  at  point  O  on  the  net  structure. 


5.  CONCLUSIONS 


Different  possibilities  for  treating  the  multiple  junction  problem  in  the  time  domain  have 
been  compared.  They  have  been  called; 

Tl.  Overlapping-segments  treatment. 

T2.  Half-segments  treatment. 

T3.  Sum  treatment. 

The  three  junction  treatments  proposed  were  specifically  applied  to  the  study  of  the 
interaction  of  a  gaussian  pulse  with: 

1. -  A  cross-wire  structure. 

2. -  A  junction  of  two  wires  with  different  radii. 

3. -  A  net  of  twelve  interconnected  wires. 

4. -  A  dodecahedron  structure. 

To  analyze  the  convergence  of  the  method  with  the  number  of  segments,  the  root  mean 
square  current  and  the  root  mean  square  error  of  the  current  induced  on  a  point  on  the  structure 
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CURRENT(mA).  IMAGINARY  PART  CURRENT(mA).  REAL  PART 


To  ascertain  their  accuracy  we  applied  the  Fourier  Transform  to  our  time  domain  results 
and  compared  them  with  some  well  known  codes  that  work  in  the  frequency  domain  (NEC  and 
AWAS). 

In  the  light  of  our  criteria,  the  overlapping  procedure  seemed  to  be  that  which  provided 
the  most  accurate  solution,  requiring  the  structure  to  be  divided  into  the  least  number  of 
segments.  This  treatment  was,  therefore,  the  one  finally  included  in  the  computer  program  we 
have  developed  for  studying  the  interaction  of  EMP  with  thin-wire  structures,  the  DOTIGl  code. 


Fig.  16  Geometry  of  the  dodecahedron  structure.  1=  0.5m.  Wire  radii  =  0.008m. 
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Fig.  17  Current  distribution  versus  time  at  point  O  on  the  dodecahedron. 


Fig.  18.  Current  magnitude  at  point  O  on  the  dodecahedron  structure  versus  frequency. 
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Optimal  Location  for  Matching  Points 
for  Wire  Modelling  with  MMP 


P.  Leuchtmann 

ETH  Zentrum,  CH-8092  Zurich,  Switzerland 

Abstract:  It  is  a  basic  premise  in  electromagnetic  field  modelling  that  the  tangential 
electromagnetic  field  on  the  surface  of  ideal  conductors  must  vanish.  When  point  matching 
is  used  to  impose  this  condition,  the  location  of  the  matching  points  must  he  chosen.  This 
paper  treats  the  special  case  of  thin  wires.  It  is  shown  that  for  common  approximations  of 
the  wire  currents,  the  optimum  locations  of  the  matching  points  are  well  defined  and  that 
forcing  the  boundary  condition  beside  these  well  defined  matching  points  would  increase 
the  errors.  For  both  piecewise-linear  and  staircase  approximation  of  the  current,  explicit 
formulae  for  the  optimum  location  of  the  matching  points  are  given. 


Introduction 

To  find  the  current  in  a  wire  one  uses  an  equation  of  the  form 

E,  =  0,  (1) 

where  Ez  is  the  longitudinal  tangential  component  of  the  electric  field  on  the  surface  of 
the  wire.  If  the  wire  is  thin,  Ez  is  small  compared  to  other  field  components  at  the 
same  location.  Nevertheless  it  is  the  only  component  to  be  matched.  For  an  approximate 
solution,  Ez  will  never  be  zero  on  the  whole  surface  of  the  wire,  but  only  in  particular 
points.  The  most  simple  procedure  for  finding  an  approximate  solution  of  (1)  is  point¬ 
matching,  where  the  wire  is  segmented  and  then  one  matching  point  is  placed  in  the 
center  of  each  wire  segment.  MMP  [1]  uses  an  extended  point  matching  technique  and  lets 
it  up  to  the  user  to  place  any  number  of  matching  points.  (A  minimum  of  one  matching 
point  per  segment  is  necessary.)  Equations  (1)  are  then  solved  in  the  sense  of  least  squares. 

It  has  been  found  that  the  results  obtained  by  the  MMP-program  are  unstable  in 
certain  cases.  In  particular,  an  increase  of  the  number  of  matching  points  did  not  lead  to 
better  results  [2].  A  closer  look  at  the  error  of  (1)  along  wires  shows  a  great  regularity  of 
the  error  distribution. 

In  this  paper,  we  shall  show  that  the  behavior  of  the  error  of  Ez  along  the  wire  is 
well  defined.  Depending  on  the  expansion  functions  of  the  current,  there  are  zeros  of  Ez 
at  certain  points.  Besides  these  points,  Ez  should  never  be  forced  to  be  zero,  since  such  a 
forcing  would  affect  the  current  in  a  highly  unintended  manner. 

The  following  scheme  gives  the  principal  way  to  find  the  results: 

•  Calculate  of  Ez  of  an  arbitrary  current  I(z)  on  a  wire  (analytically) 

•  Approximate  I(z)  by  a  staircase  function  (like  MININEC)  or  a  piecewise-linear  func¬ 
tion  (MMP)  and  calculate  Ez  of  the  approximated  current  (analytically) 

•  Define  the  systematic  error  as  the  difference  —  Ez 

•  Use  zeros  of  the  systematic  error  as  optimum  locations  for  matching  points 
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The  reason  why  the  zeros  of  the  systematic  error  are  optimum  locations  for  matching 
points  is  quite  obvious:  if  one  would  try  to  push  down  the  error  at  a  location,  where  the 
systematic  error  is  nonzero,  the  error  of  the  solution  is  increased,  since  the  systematic  error 
ccLnnot  be  eliminated. 

As  a  principal  result,  we  shall  find  that  these  optimum  locations  are  highly  independent 
of  the  global  shape  of  I(z),  but  only  a  consequence  of  the  local  ‘misfit’  of  the  expansion 
functions  vs.  the  exact  field.  As  a  matter  of  fact,  the  optimum  locations  of  the  matching 
points  depend  on  both  the  expansion  functions  and  the  unknown  field.  Therefore,  both  of 
them  have  to  be  considered  to  find  good  rules  for  the  placing  of  the  matching  points.  To 
get  rid  of  this  highly  contradictory  situation  —  the  unknown  field  must  be  known  —  we 
shall  consider  that  class  of  unknown  current  distributions  which  can  be  represented  by  a 
converging  Taylor  series.  It  is  worth  noting  that  the  behavior  of  the  expansion  functions 
alone  is  not  sufiicient  to  find  the  optimum  locations  of  the  matching  points. 


Taylor-Representation  of  the  Exact  Current-  and  Charge  Distribution 

Let  us  consider  a  straight  (thin)  ^r-oriented  wire  with  radius  r  which  is  carrying  a  current 
Iexa.ct{z).  On  a  piece  of  finite  length  L  of  this  wire,  the  current  may  be  represented  by  a 
power  series^ 

N 

•ifexact(^)  «  ^  (2) 

n=0 

We  assume,  that  the  finite  power  series  (2)  is  able  to  represent  the  function  /exact(^)  within 
the  length  L  so  that  no  further  error  considerations  are  necessary.  Later  on,  we  shall  divide 
up  the  length  L  into  segments  and  we  shall  focus  on  the  local  effects  of  this  segmentation. 
From  this  point  of  view,  the  assumptions  are  reasonable. 

Due  to  the  charge  conservation  principle,  the  current  distribution  I(z)  is  accompanied 
by  the  charge  distribution^ 


iuj  dz 


n=l 


nanZ^  ^ 


(3) 


As  a  further  assumption,  we  state,  that  L  is  much  smaller  than  the  wavelength  A,  so 
that  the  electric  field  E  may  be  calculated  from  q{z')  by  means  of  electrostatics.  Since  we 
are  only  interested  in  local  effects,  we  do  not  consider  the  field,  which  is  produced  by  the 
currents  and  charges  beyond  the  length  L.  This  field  is  not  necessarily  small  but  almost 
constant  along  L.  In  particular,  it  does  not  vary  on  a  segment  of  L.  For  this  reason,  it 
has  no  importance  for  the  location  of  the  matching  points  relative  to  the  segment. 

The  ^-component  of  the  S-field  is  obtained  from  q{z)  by  Coulomb’s  integral: 


Ezo{r,z)  = 


Vr*  +  (z  -  z'f 


L 


{zr-\z-z') 

i/r2  +  {z-  z'y^ 


(4) 


^  Convergence  is  given,  if  L  is  sufficiently  small  and  I{z)  is  smooth. 
^  Time  dependency  of  is  assumed! 
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Figure  1  Diiferent  current  distributions  cause  different  distributions  of  the  z-component 
of  E,  The  pictures  show  the  four  parts  of  Ez  for  n  =  1 . .  .4  with  otn  =  1, 
r  =  0,1  and  L  =  10.  At  the  middle  and  at  the  bottom,  details  of  the  interval 
0  ^  z  <  2  are  shown.  The  numerical  values  are  obtained  from  (3),  omitting  the 
factor  l/(47r«a;£:).  We  focus  on  the  area  around  z  =  0. 


The  integral  on  the  right  hand  side  may  be  solved  analytically  for  any  n.  These  solutions 
and  the  plots  of  these  functions  have  been  obtained  using  Stephen  Wolfram’s  ‘mathemat- 
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ica’.  [3]  (See  fig.  l!) 

Due  to  the  superposition  principle,  each  term  may  be  treated  alone.  Provided  that 
the  optimum  locations  for  matching  points  being  the  same  for  each  n,  these  locations  are 
optimum  for  any  I{z).  This  is  exactly  what  we  shall  find. 

At  the  ends  of  we  assume  a  constant  continuation  of  the  current.  Otherwise,  the 
steps  in  the  current  distributions  would  cause  point  charges  at  the  ends,  which  would 
disturb  the  results  in  producing  an  additional  large  electric  field. 


Staircase  Approximation 

The  wire  of  length  L  is  cut  into  J  segments  of  equal  length  I  =  L/  J.  Then  the  function 
I(z)  is  approximated  by  a  staircase.  (See  fig.  2!)  The  same  approximation  is  used  in 
MININEC  [4].  Note  that  the  point-matching  technique  which  is  described  in  this  paper  is 
different  to  the  (more  sophisticated)  procedure  used  in  MININEC. 


piecewise-Iinear  functions  (MMP).  Each  type  of  approximation  is  accompa¬ 
nied  by  its  typical  (erroneous)  charge  distribution.  It  is  this  charge  distribution 
which  causes  a  relatively  big  error  in  the  local  distribution  of  the  longitudinal 
electric  field. 

Let^ 

(5) 

be  the  z-coordinate  of  the  boundary  between  the  j-th  and  the  {j  -f-  l)-th  segment.  In  the 


3 


For  simplicity’s  sake,  J  is  assumed  to  be  even! 
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j-th  segment  the  current  is  constant  and  equal  to 


for  -  f  +  1  <  j  <  ^ 

ioTj  =  -j  (6) 

U(^j-i)  fori  =  1  +  1 

I.e.,  we  assume  the  approximation  being  exact  in  the  center  of  the  segment.  The  special 
definition  of  the  last  two  currents  is  due  to  the  constant  continuation  beyond  L.  This 
current  distribution  leads  to  point  charges  Qj  located  at  the  segment  boundaries  Zj.  One 
has 

Qi  =  -  Ij)  (7) 

From  that,  we  obtain  for  the  ^-component  of  the  resulting  electric  field 


Ezsir,  z) 


1  ^ 
=  J-  y- 

A—/ 


Qj  •  -  ^j) 

/r2  +  (z  -  Zj} 


{Ij+l  -  Ij){z  -  Zj) 


4:TriLje 


j=-i 


The  sum  has  J  +  1  terms,  because  a  line  of  J  segments  has  J  +  1  boundaries.  Figure  3 
shows  the  big  variation  of  along  the  wire  and  in  particular  along  a  single  segment. 

Let  us  define  now  the  systematic  error  as  the  difference  between  the  “exact”  result 
(4)  and  the  approximation  (8): 


F:(r,z)  =  E^,(r,z)-EUr,z)  (9) 

where  the  upper  index  n  denotes  the  power  of  the  current:  I{z)  =  2;”.  It  is  clear,  that 
Fg  must  vanish,  since  the  staircase  approximation  of  a  constant  function  is  exact.  Figure 
4  shows  the  functions  F”  along  the  segments  1  and  2.  The  values  around  z  =  0  are 
not  typical,  because  of  the  superimposition  of  a  non  local  effect  due  to  the  symmetrical 
situation:  With  odd  powers  of  the  current  distribution  (i.e.,  even  powers  for  the  charge 
distribution),  the  far  action  of  the  discretization  errors  are  cancelled  due  to  symmetry. 
However,  with  even  powers  of  the  current  distribution  (i.e.,  odd  powers  for  the  charge 
distribution),  these  actions  add  together  and  cause  a  vertical  shift  of  the  curve.  For  these 
reasons,  the  behavior  in  the  second  segment  is  more  typical. 

In  general  we  can  say  that  for  all  powers  n,  the  systematic  error  has  a  relatively  fiat 
zero  in  the  middle  of  the  segment.  For  that  reason,  it  is  optimum  to  place  the  matching 
points  in  the  center  of  the  segment.  The  second  zero  (close  to  the  segment  boundaries) 
is  not  useful  for  this  purpose,  because  the  error  function  is  much  too  steep  in  this  area. 
Considering  the  exact  values  of  Ez  (see  fig.  1!),  we  find  that  the  error  is  not  negligible  at 
all. 

As  a  preliminary  conclusion  we  state:  For  staircase  approximation  it  is  optimum  to 
place  matching  points  in  the  center  of  the  segments  —  and  nowhere  else.  This  is  due  to 
the  fact  that  the  systematic  error  cannot  be  eliminated.  A  placing  of  a  matching  point  at 
a  position  with  a  high  systematic  error  would  increase  the  global  error  of  the  field  in  such 
a  sense,  that  the  sum  of  the  two  vanishes. 
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Figure  3  The  staircase  approximation  of  the  current  gives  a  much  less  smooth  behavior 
of  Ez  than  the  smooth  current  in  fig,  1.  See  particularly  the  big  values  in  the 
center  (near  z  =  0)!  The  picture  shows  the  four  parts  of  Ez  for  n  =  1 . ,  .4  with 
an  =  1,  r  =  0.1,  L  =  10  and  J  =  10  (— >  segment  length  I  =  1).  At  the  middle 
and  at  the  bottom,  details  of  the  interval  0  <  z  <  2  are  shown.  The  numerical 
values  are  obtained  from  (9),  omitting  the  factor  l/(47rza;6:). 
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Figure  4  The  systematic  error  Fg  of  Ez  (using  a  staircase  approximation  of  the  current) 
depends  on  the  wire  radius  r.  Since  the  length  of  the  segments  has  been  kept 
constant,  the  ratio  length/radius  of  a  segment  is  actually  varied.  The  change  of 
the  vertical  scale  at  the  bottom  of  each  picture  shows  the  poor  flatness  of  the 
zero.  As  expected,  the  absolute  value  of  the  errors  decreases  with  increasing 
radius.  Note  that  the  numerical  values  in  fig.  1  are  related  to  a  wire  radius  of 
r  ~  0.1. 


27 


Figure  A  Continued,  see  above! 


Piecewise-linear  Approximation  (MMP) 

As  it  has  been  already  done  in  the  previous  section,  the  wire  of  length  L  is  split  up 
into  J  segments  of  equal  length  I  =  Lj  J .  Other  than  before,  the  function  I(z)  is  now 
approximated  by  a  piecewise-linear  function*^.  In  the  j-th.  segment,  the  current  is  a  linear 
function: 

In  the  MMP-programs,  I{z)  is  approximated  using  the  functions  coskz  and  sin^0,  where  k  is  the 
wave  number  of  the  surrounding  medium  [2].  Since  the  segments  are  small  compared  to  wavelength,  these 
two  functions  tend  to  1  and  z. 
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i.e.,  we  assume  that  the  approximation  is  exact  at  the  segment  boundaries.  (See  fig.  2!) 
This  approximation  of  the  current  distribution  leads  to  a  constant  charge  distribution  qj 
on  the  j-th  segment: 


<lj  = 


1 


liO  I 

and  the  ^-component  of  the  resulting  electric  field  is 


(11) 


Ezi{r,z) 
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Obviously,  this  sum  has  one  term  less  than  the  one  in  (8).  Figure  5  shows  the  variation  of 
Ez  which  is  now  much  smaller  than  the  one  shown  in  fig.  3. 

Again,  the  systematic  error  F”  is  defined  as  the  difference  between  the  “exact”  result 
(4)  and  the  approximation  (12): 


Fr{r,z)  =  E^,(r,z)-E^„(r,z)  (13) 

where  the  meaning  of  the  upper  index  is  the  same  as  in  (9).  It  is  clear,  that  both  F°  and 
Fl  vanish,  since  the  piecewise-linear  approximation  of  a  linear  function  is  exact. 

Considering  the  systematic  error  we  find  —  as  in  the  staircase  approximation  —  a 
great  regularity,  which  is  qualitatively  equal  for  all  powers.  On  the  other  hand,  there  is 
a  fundamental  difference:  The  zeros  of  the  systematic  error  do  not  occur  in  the  middle  of 
the  segments,  but  more  close  to  the  segment’s  boundaries.  (See  figs.  6  and  7!) 

Due  to  the  global  effect  in  the  middle  of  L  (vertical  shift  of  the  even  order  curves  and 
a  zero  at  z  =  0  for  the  odd  orders  due  to  symmetry),  the  first  zero  after  z  =  0  is  not 
used.  (For  the  numbering  of  the  zeros,  see  fig.  6  at  the  bottom.)  The  exact  location  of  the 
zeros  of  the  systematic  error  is  a  function  of  the  radius  r,  the  order  n  and  the  number  of 
the  zero.  Figure  8  shows  this  locations  as  a  function  of  the  wire  radius  r.  Note  the  great 
regularity! 

Due  to  the  fact  that  the  lowest  order  (here  n  =  2)  is  most  important  and,  on  the 
other  hand,  curves  for  different  zero  numbers  axe  —  at  least  for  small  radii,  i.e.  thin  (long) 
segments  —  almost  equal,  it  makes  sense  to  use  e.g.  the  third  zero  of  the  second  order  as 
a  reference.  The  bigger  variation  at  r //-ratios  more  than  1/3  does  not  matter,  since  the 
absolute  values  of  the  systematic  errors  decrease  with  increasing  wire  radius. 

We  conclude:  For  piecewise-linear  approximation,  matching  points  in  the  center  of  the 
segments  should  be  avoided.  There  are  two  optimal  locations  for  each  segment,  symmet¬ 
rically  placed  at  a  distance  d  away  from  the  segment’s  boundaries. 
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gives  an  unquiet  behavior  of  Ez^  but  the  improvement  with  respect  to  the  stair¬ 
case  approximation  (fig.  3)  is  obvious.  The  pictures  show  the  four  parts  of  Ez 
for  n  =  1 . .  .4  with  an  =  1,  r  =  0.1,  L  =  10  and  J  =  10  (-^  segment  length 
I  =  1).  At  the  middle  and  at  the  bottom,  details  of  the  interval  0  <  z  <  2 
are  shown.  The  numerical  values  are  obtained  from  (12),  omitting  the  factor 
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Figure  6  The  systematic  errors  Ff'  (in  the  case  of  a  piecewise-linear  approximation  for 
the  current)  are  very  similar  for  different  powers.  At  the  top,  the  figure  shows 
the  systematic  errors  of  the  order  n  =  2..  .4.  At  the  bottom,  the  most  im¬ 
portant  order  2  is  given  individually.  The  numbering  of  the  zeros  (numbers  in 
circles)  starts  at  the  second  one.  This  is  due  to  the  symmetry  caused  fact,  that 
the  even  curves  (from  4-th  order  on)  are  vertically  shifted  down  around  z  =  0, 
while  odd  curves  have  a  zero  at  z  =  0.  The  distance  between  the  zeros  and  the 
segment  boundaries  (here  integer  valued  z’s!)  is  the  important  parameter  for 
placing  matching  points. 
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of  the  current  depends  on  the  wire  radius  r.  Since  the  length  of  the  segments 
has  been  kept  constant,  the  ratio  length/radius  of  a  segment  is  actually  varied. 
Note  that  the  numerical  values  in  fig.  1  are  related  to  a  wire  radius  of  r  —  0.1. 
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A 

Figure  7  Continued,  see  above! 
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Polynomial  Approximation  of  the  Zero’s  Relative  Position 

We  have  chosen  as  a  reference  the  zero  number  3  of  the  quadratic  term  of  the  current.  It  is 
located  at  a  distance  dz  away  from  the  segment’s  boundary.  This  distance  dz  is  a  function 
of  the  wire  radius  r:  dz  =  dz(r).  ‘Mathematica’  [3]  does  the  job  to  expand  the  function 
dz(r)  into  a  power  series  of  m-th  order.  We  normalize  the  segment  length  to  1  and  find: 

dl{r)  =  0.1246832  +  0.572177r  -  0.7090428r2 

df  (r)  =  0.1085357  +  0.893254r  -  2.043528r2  +  1.480294r^  (14) 

dj{r)  =  0.0978577  +  1.245951r  -  4.679222r2  +  8.298414r^  -  5.672312r^ 
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The  pictures  show  the  distance  di  between  the  i-th  zero  and  the  next  segment 
boundary  as  a  function  of  the  wire  radius  r.  Only  the  second  and  the  third  or¬ 
der  is  shown.  The  position  of  the  zeros  with  higher  numbers  are  inaccurate  due 
to  end  effects  of  the  model  (only  10  segments). 

The  upper  index  denotes  the  order  of  the  approximations,  which  are  valid  within  the 
interval  0.001  <  r  <  0.6. 

The  deviation  of  the  true  value  is  maximum  at  the  border  of  this  interval.  We  denote 
the  biggest  deviation  downwards  (approximation  smaller  than  true  value)  with  Im  and  the 
biggest  deviation  upwards  (approximation  bigger  than  true  value)  with  Um-  Mathematica 
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delivers 


h  =  0.05406: 
h  =  0.03824: 
h  =  0.02791: 
h  =  0.02208: 


U2  =  0.00926 
U3  =  0.00634 
U4  =  0.00471 
Us  =  0.00396 


(15) 


Comparing  these  numbers  with  the  true  value  (?»  0.2)  and  considering  the  variation  be¬ 
tween  the  different  zeros  (see  fig.  8!),  an  approximation  of  ds^r)  of  4-th  order  is  sufficient. 


Summary 

We  have  found  that  the  positions  of  matching  points  along  thin  wires  must  be  chosen 
carefully.  Depending  on  the  segmentation  and  the  type  of  the  approximation  of  the  currents 
in  the  wire,  optimum  matching  point  positions  are  different.  The  results  have  been  found 
by  examining  the  error  produced  by  a  ‘best  fit’  of  some  expansion  set  (staircase  and 
piecewise-linear  expansion  functions)  to  an  arbitrary  current  distribution. 

The  free  choice  of  matching  points  in  the  MMP-code  is  an  error  source  which  could  be 
avoided.  Using  MMP  one  should  set  two  matching  points  for  each  segment,  symmetrically 
located  at  a  distance  d  away  from  the  segment’s  boundaries  (see  eq.  (14)!).  This  choice  is 
optimum  at  all  locations  where  the  current  I{z)  is  either  concave,  convex  or  linear. 

As  a  further  result,  we  find  that  the  point  matching  technique  is  preferred  to  any  other 
method  to  satisfy  boundary  condition  (1),  as  long  as  the  current  is  a)  approximated  by  fixed 
functions  and  b)  strictly  coupled  to  the  electric  field  through  Maxwell’s  equations.  This 
second  point  is  not  fulfilled,  e.g.,  by  MININEC:  it  uses  a  piecewise-constant  approximation 
of  the  charge  distribution,  which  would  be  related  to  a  piecewise-linear  current  distribution. 
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ABSTRACT 


This  paper  deals  with  the  numerical  modelling  of  the  surface  of  a  structure  with  a 
wire  grid  or  a  surface  patch  for  Electromagnetic  (EM)  interaction  analysis.  Surface 
currents  and  fields  on  a  wire  grid  model  are  computed  using  the  Numerical 
Electromagnetic  Code  (NEC).  The  results  are  compared  with  those  obtained  on  a 
triangular  surface  patch  model  using  an  Electric  Field  Integral  Equation  (EFIE)  formulation. 
Simple  structures  such  as  a  square  plate  as  well  as  complicated  structures  such  as  an 
aircraft  are  considered.  Good  agreement  is  obtained  in  most  cases. 


INTRODUCTION 

Numerical  modelling  of  a  surface  with  a  wire  grid  has  been  used  for  many 
electromagnetic  antenna  radiation  and  scattering  problems  [1-4].  It  has  many  attractive 
features  and  is  capable  of  giving  reliable  far-field  results  if  a  proper  choice  of  wire 
diameters  and  grid  size  is  made  [5].  However,  some  questions  have  been  raised 
regarding  the  validity  of  using  wire  grid  models  of  a  closed  surface  when  calculating 
electromagnetic  pulse  (EMP)  interactions  [6].  These  questions  arise  because,  for 
calculating  EMP  interaction,  one  must  determine  the  currents  and  current  densities 
induced  on  the  surface.  In  the  case  of  a  structure  with  a  closed  surface,  the  field  must 
vanish  identically  in  the  interior  region,  whereas  for  the  wire  grid  model,  there  is  an 
evanescent  reactive  field  clinging  to  both  sides  of  the  grid.  Even  for  a  structure  with  an 
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open  surface  (e.g.  a  plate)  there  are  questions  regarding  the  variation  of  the  current 
around  the  wire  circumference,  and  regarding  the  behaviour  of  the  fields  between  the 
wires  and  at  the  junctions. 

There  are  also  questions  regarding  the  choice  of  the  wire  diameters  in  a  wire  grid 
representation.  There  are  rules  for  choosing  wire  diameters  for  obtaining  reliable  results 
for  far  field  quantities.  It  is  not  clear  whether  the  same  rules  apply  for  the  estimation  of 
currents  and  current  densities  induced  on  a  structure  exposed  to  EMP. 

This  paper  attempts  to  answer  some  of  these  questions  by  comparing  wire  grid 
and  surface  patch  modelling  for  a  number  of  structures.  Simple  structures  such  as  a 
square  plate  as  well  as  complicated  structures  such  as  an  aircraft  are  considered.  The 
Numerical  Electromagnetic  Code  (NEC)  [1]  is  used  for  the  wire  grid  models  and  the 
Electric  Field  Integral  Equation  (EFIE)  [7]  is  used  for  the  surface  patch  models.  Surface 
currents  and  current  densities  induced  on  a  structure  due  to  an  incident  electromagnetic 
wave  are  computed  using  the  two  techniques  and  compared.  Some  modifications  are 
required  in  the  existing  codes  in  order  to  deal  with  large  number  of  wire  or  patch 
segments  and  to  overcome  difficulties  involved  in  modelling  surfaces  where  three  patches 
have  one  common  edge.  Although  the  comparison  is  done  in  the  frequency  domain,  the 
time  domain  behaviour  can  be  obtained  using  the  Fourier  transform. 


PROCEDURE 

A  wire  grid  representation  of  the  surfaces  of  the  body  which  is  to  be  analyzed  is 
first  made  following  the  rules  given  in  the  NEC  documentation  [1].  For  a  simple  structure 
such  as  a  square  plate  this  is  fairly  straightforward.  However,  for  a  complicated  structure 
the  creation  of  a  wire  grid  model  is  difficult  and  one  must  use  specialized  software  such 
as  DIDEC.DREO  [8].  The  diameter  of  the  wires  is  chosen  to  follow  the  criterion  that  the 
surface  area  of  the  wire  parallel  to  one  linear  polarization  should  be  equal  to  the  surface 
area  of  the  solid  surface  being  modelled  [9].  The  NEC  is  then  used  to  determine  currents 
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at  the  centre  of  each  wire  segment.  The  induced  H-field  or  surface  current  density  at  a 
point  on  the  equivalent  surface  represented  by  the  wire  is  then  computed  by 

^induced  “  l/27r3, 

where  I  is  the  current  induced  in  the  wire  obtained  by  NEC  computations  and  a  is  the 
radius  of  the  wire. 

In  order  to  check  the  accuracy  of  the  wire  grid  representation  for  determining 
surface  current  and  currents,  a  surface  patch  model  is  created.  This  surface  patch  model 
uses  the  Electric  Field  Integral  Equation  (EFIE)  instead  of  the  magnetic  field  integral 
equation  (MFIE)  as  proposed  in  the  NEC.  The  surface  of  the  body  is  modelled  by 
triangular  patches  and  a  modified  version  of  the  program  EFIE  [10]  is  used  to  compute 
the  surface  current  distribution.  The  modifications  to  the  EFIE  include  the  use  of  a  new 
matrix  LU  factorization  routine  and  the  addition  of  a  solution  polisher.  The  new  matrix 
factorization  routine  results  in  a  great  reduction  in  time.  The  solution  polisher  gives  some 
extra  significant  digits  whenever  the  impedance  matrix  is  ill-conditioned.  The  induced  H- 
field  output  of  the  EFIE  is  available  at  the  edges  of  every  triangular  patch  and  is  directed 
perpendicular  to  the  edge  and  the  normal  to  the  surface. 

Since  the  NEC  wire  grid  modelling  computations  yield  current  and  induced  H-field 
along  the  axis  of  a  wire  and  the  EFIE  triangular  surface  patch  modelling  computations 
yield  induced  H-field  perpendicular  to  the  triangle  edge,  care  has  to  be  taken  in 
comparing  data  from  the  two  models.  The  need  for  this  care  is  demonstrated  in  Figures 
1a  to  1c  which  show  the  wire  grid  and  the  equivalent  surface  patch  models  for  a  square 
plate.  The  Figures  show  90  wire  rectangles  for  the  wire  grid  representation  and  84 
triangular  patches  for  the  surface  patch  representation.  This  discretization  is  chosen  to 
allow  easy  comparison  of  induced  fields  or  charge  densities  along  two  principal  cut  lines 
of  the  square  plate.  Other  discretizations  are  permitted  as  long  as  the  rules  of  EFIE  and 
NEC  are  not  violated.  The  points  along  the  two  axes  (XX’  and  YY’)  at  which  the  surface 
fields  are  evaluated  on  the  wire  grid  and  the  surface  patch  representation  are  shown  in 
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the  Figures.  These  points  are  not  identical  but  can  be  made  so  if  required  without  much 
difficulty.  It  is  obvious  that  making  a  choice  of  points  for  comparison  also  presents  little 
difficulty  in  this  simple  case.  However  the  choice  of  points  for  making  a  proper 
comparison  between  two  representations  is  not  always  straightforward  for  a  complicated 
structure  such  as  an  aircraft. 

Figure  2a  shows  the  current  density  distribution  along  two  principal  cut  lines  (XX’ 
and  YY’  in  Figure  1)  on  a  1.0  wavelength  square  plate,  illuminated  by  a  plane 
electromagnetic  wave  incident  normally  on  the  plate.  In  Figure  2a  and  the  subsequent 
Figures  the  induced  current  density  normalized  by  the  incident  magnetic  field  intensity  is 
plotted  on  the  vertical  axis.  Because  of  the  normalization  this  quantity  has  no  units.  The 
comparison  between  the  wire  grid  and  triangular  surface  patch  results  is  very  good 
except  at  the  edges  of  the  plate.  The  normalized  current  densities  obtained  by  the  wire 
grid  model  are  smaller  at  the  edges  because  the  equivalent  surface  area  is  larger  and  the 
"equal-area"  rule  can  no  longer  be  followed. 

If  the  diameter  of  the  wires  at  the  edges  is  reduced  by  half  (Figure  1c),  the  wire 
grid  results  seem  to  match  very  well  with  the  surface  patch  results  even  at  the  edges. 
This  is  shown  in  Figure  2b.  The  wire  grid  results  at  the  edges  of  the  plate  perpendicular 
to  the  y-axis  are,  however,  still  not  correct.  Because  of  discontinuity,  the  normalized 
current  density  function  is  singular  (i.e.  has  an  infinite  value)  at  these  edges  [7].  This 
infinite  value  can  only  be  approached  by  the  wire  grid  results  in  the  limit,  when  the  edge 
wire  diameter  and  the  area  it  represents  become  infinitesimally  small.  The  surface  patch 
results  are  not  obtained  at  these  edges. 

Figure  3b  shows  the  dominant  current  density  distribution  along  a  cut  (line  XX’) 
through  the  symmetry  plane  of  a  1.0  wavelength  bent  square  plate  (Figure  3a).  The 
distance  S  shown  in  Figures  3b  and  3c  is  measured  along  the  plate  from  point  X  indicated 
in  Figure  3a.  The  bend  is  located  at  two-third  the  plate  width  from  this  point  and  is  parallel 
to  the  plate  edge.  A  plane- wave  with  the  E-field  parallel  to  the  bend  is  incident  normal  to 
the  large  section  of  the  plate.  The  smaller  section  is  bent  at  an  angle  of  50°  towards  the 
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shadow  side  of  the  plate.  There  are  72  (or  6x6x2)  triangular  patches  for  the  EFIE  solution 
versus  156  (or  13x12)  wire  rectangles  for  the  wire  gird  model  NEC  solution.  The  diameter 
of  the  wires  is  chosen  to  satisfy  the  equal-area  rule  discussed  earlier  (although  for 
simplicity  they  are  not  shown  in  Figure  3a).  The  comparison  of  the  wire  grid  and  the 
triangular  surface  patch  models  is  again  very  good  except  close  to  the  edges  of  the  plate. 

As  discussed  before,  the  discrepancy  results  from  the  non-compliance  of  the 
equal-area  rule  at  the  edges.  Again,  if  the  diameter  of  the  wires  at  the  edges  is  reduced 
by  half,  the  agreement  between  the  wire  grid  and  the  surface  patch  results  is  improved. 
This  is  shown  in  Figure  3c.  Earlier  comments  on  the  singularity  of  the  induced  current 
density  at  the  edges  also  apply  here. 

Figure  4a  shows  an  example  of  a  complicated  wire  grid  model  of  an  aircraft  that 
is  analyzed  here.  Figure  4b  shows  the  equivalent  surface  patch  model.  The  wire  grid 
structure  has  already  been  analyzed  for  scattering  and  radiation  properties  [9].  In  our 
case,  the  wire  grid  structure  consists  of  326  wire  segments  and  the  surface  patch 
representation  consists  of  398  triangular  patches.  A  plane  wave  with  the  electric  field 
polarized  parallel  to  the  long  axis  of  the  aircraft  is  incident  normal  to  the  top  of  the  aircraft. 
The  x-axis  runs  parallel  to  the  length  of  the  aircraft  (hence  parallel  to  the  incident  electric 
field),  with  the  reference  point  x=0  located  in  the  middle  of  the  attachment  of  the  wings 
to  the  fuselage,  and  with  positive  x-coordinate  values  towards  the  front  of  the  aircraft. 
Figure  5a  compares  the  normalized  current  densities  of  the  two  models  along  the  top 
fuselage  of  the  aircraft  at  20  MHz.  The  comparison  shows  a  good  agreement  between 
the  surface  patch  and  the  wire  grid  models.  Figure  5b  compares  the  normalized  current 
density  around  the  circumference  of  the  fuselage  at  the  front  of  the  aircraft  at  10  MHz. 
The  agreement  between  the  wire  grid  and  the  surface  patch  models  is  again  quite  good. 

The  discrepancy  is  partly  explained  by  the  fact  that  the  points  of  observation  are 
not  identical  in  this  case.  In  a  complicated  structure  such  as  an  aircraft,  it  is  difficult  to 
maintain  the  same  observation  points  for  the  wire  grid  and  the  surface  patch  models. 
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Figure  6  shows  the  wire  grid  and  the  triangular  surface  patch  models  of  a  metallic 
cube  with  an  open  top.  Each  side  of  the  cube  is  1.0  wavelength  long.  The  models  consist 
of  405  wire  rectangles  and  408  triangular  surface  patches,  respectively.  A  polarized  plane 
wave  is  incident  on  the  open  face  as  shown.  Figure  7  shows  the  normalized  induced 
current  density  along  the  lines  XX’  and  ZZ’  on  one  face  of  the  cube  for  the  two  models. 
The  agreement  between  the  models  is  very  good  even  at  the  edges  of  the  cube.  It  is  to 
be  noted  that  the  wire  diameters  at  the  top  edges  (only  at  the  open  end)  parallel  to  the 
incident  field  have  been  reduced  to  follow  the  equal-area  rule.  The  current  density  at  the 
open  edge  again  exhibits  a  singular  behaviour,  and  the  earlier  comments  on  the  wire  grid 
results  apply. 


EMP  INTERACTION  ANALYSIS 

The  above  analysis  can  be  performed  for  a  number  of  frequencies  as  long  as  the 
rules  of  the  wire  grid  and  those  of  the  surface  patch  representation  are  followed. 
Response  to  a  pulse  excitation  such  as  the  EMP  can  be  obtained  by  convolution  with  the 
EMP  spectrum  and  taking  a  fast  Fourier  transform  of  the  response  at  a  proper  number 
of  frequencies  [11].  The  procedure  is  simple  for  the  cases  of  the  flat  square  plate,  the 
bent  square  plate,  and  the  open-topped  cube.  However,  for  the  case  of  the  aircraft, 
difficulties  arise  because  of  the  large  number  of  wire  segments  or  triangular  patches 
required  at  the  upper  end  of  the  EMP  spectrum.  The  CPU  time  required  in  this  case  is 
considerable  even  for  a  single  frequency. 


CONCLUSIONS 


We  have  computed  induced  surface  fields  for  a  number  of  structures  for  plane 
wave  incidence.  The  limited  number  of  examples  treated  indicate  that  the  results  obtained 
by  the  wire  grid  and  the  triangular  surface  patch  models  show  good  agreement  and  both 
models  may  be  used  to  compute  induced  surface  fields. 
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Figure  1a.  Wire  grid  (NEC)  model  of  a  flat  square  plate  with  equal- 
area  wire  diameters. 


Figure  1b.  Surface  patch  (EFIE)  model  of  the  flat  square  plate, 
equivalent  to  the  wire  grid  model  of  Fig.  la. 


46 


Figure  1c.  Wire  grid  (NEC)  model  of  the  flat  square  plate  with  edge 
wire  diameters  halved. 


SQUARE  FLAT  PLATE 


Figure  2a.  Distribution  of  dominant  component  of  current  density 
along  XX’  and  YY’  on  1 .0  lambda  plate;  x  and  solid  iine  =  wire  grid 
(NEC)  model  of  Fig.  la,  o  =  surface  patch  (EFIE)  model  of  Fig.  1b. 
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Figure  3a.  Surface  patch  and  wire  grid  models  of  a  bent  square 
plate 
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BENT  PLATE 
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Figure  3b.  Distribution  of  dominant  component  of  current  density 
along  XX’  on  the  bent  square  plate  of  Fig.  3a;  x  and  solid  line  =  wire 
grid  (with  equal-area  wire  diameters),  o  =  surface  patches. 


BENT  PLATE 
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Figure  3c.  Same  as  Fig.  3b,  but  the  wire  grid  model  has  edge  wire 
diameters  halved. 


Figure  4b.  Equivalent  surface  patch  model  of  the  aircraft  in  Fig.  4a. 
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AIRCRAFT  TOP  FUSELAGE 


Figure  5a.  Normalized  current  density  distribution  along  the  top  of 
the  fuselage  of  the  aircraft  in  Fig.  4  at  20  MHz;  solid  line  =  wire  grid 
(NEC),  dashed  line  =  equivalent  surface  patches  (EFIE). 


AIRCRAFT  AROUND  FRONT  FUSELAGE 
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Figure  5b.  Normalized  current  density  distribution  around  the  front 
fuselage  of  the  aircraft  in  Fig.  4  at  10  MHz;  solid  line  =  wire  grid 
(NEC),  dashed  line  =  surface  patches  (EFIE). 
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Figure  6.  The  wire  grid  and  surface  patch  models  of  an  open-topped 
cube. 


OPEN-TOPPED  CUBIC  BOX 


Figure  7.  Normalized  current  density  along  XX’  and  ZZ’  on  the  open- 
topped  cube  of  Fig.  6;  x  and  solid  line  =  wire  grid  (NEC),  o  = 
surface  patches  (EFIE). 
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ON  THE  FUNCTIONING  OF  A  HELICOPTER-BORNE  HF  LOOP  ANTENNA 
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ABSTRACT:  Numerical  techniques  allow  engineers  to  evaluate  the  performance  of  antennas 
on  complex  structures.  These  techniques  can  provide  valuable  physical  insights  into  the  overall 
functional  performance  of  such  antennas.  This  short  paper  reports  on  the  use  of  NEC2  to 
investigate  the  radiation  characteristics  of  a  loop  antenna  mounted  below  the  tailboom  of  a 
helicopter  in  the  frequency  band  2-15  MHz.  It  was  concluded  that  such  a  loop  antenna  serves 
mainly  to  excite  a  dominant  electrical  dipole  mode  of  operation  for  frequencies  less  than  the 
lowest  natural  electrical  resonance  frequency  of  the  helicopter  itself  but  greater  than  some 
frequency  near  the  low  end  of  the  HF  band.  A  limited  set  of  measurements  performed  on  a 
scale  model  of  a  helicopter  generally  supports  the  conclusions  drawn  from  the  numerical  results 
predicted  using  NEC2.  The  reported  result  is  of  importance  in  so-called  'nap  of  the  earth’  HF 
communications  from  helicopters. 

1.  INTRODUCTION 

The  method  of  moments  [1]  is  used  extensively  in  the  design  and  analysis  of  wire 
antennas  as  well  as  the  evaluation  of  the  performance  of  such  antennas  on  a  variety  of 
vehicles  and  structures  [2].  The  method  has  proved  to  be  extremely  useful  for  the 
evaluation  of  the  performance  of  HF  (3-30  MHz)  antennas  on  aircraft,  because  of  problems 
associated  with  measuring  the  overall  performance  either  during  flight  tests  or  using  scale 
models.  For  numerical  analysis  the  aircraft  fuselage  and  wings  are  modelled  as  wire-grid 
structures.  The  method  is  useful  for  examining  the  detailed  performance  of  aircraft 
antennas  in  the  HF  band  and  as  an  engineering  tool  in  designing  such  antennas. 

Kubina  [3]  presents  guidelines  for  numerically  modelling  wire  grid  structures.  Along  with 
other  examples  he  also  gives  predicted  and  measured  radiation  patterns  for  a  loop  antenna 
on  a  CHSS-2/Sea  King  helicopter  at  frequencies  of  2.6, 4.1,  6.0  and  8.1  MHz.  The  loop  was 
mounted  on  the  tailboom  in  the  horizontal  (x-y)  plane.  Owen  [4]  also  gives  guidelines  for 
numerical  modelling  and  presents  measured  and  predicted  radiation  patterns  for  a 
helicopter  mounted  loop  antenna  at  a  frequency  of  8  MHz.  In  this  case  the  loop  antenna 
was  mounted  in  the  vertical  (x-z)  plane  below  the  tailboom  of  the  helicopter.  In  [3]  and  [4] 
the  emphasis  appeared  to  be  on  the  validation  of  computed  results,  and  how  best  to  achieve 
good  comparisons  with  the  measured  results  by  ’tuning’  the  numerical  model.  Neither 
author  examined  the  various  radiation  modes  which  could  influence  the  radiation  patterns. 
Burberry  [5]  points  out  that  there  are  two  modes  of  radiation  present,  that  due  to  the 
antenna  itself,  as  well  as  an  additional  mode  due  to  the  currents  induced  on  the  airframe 
by  coupling  with  the  antenna.  From  an  examination  of  the  results  presented  in  [3]  and  [4] 
this  author  concluded  that  a  longitudinal  electrical  dipole  mode,  corresponding  with  the 
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longest  dimension  of  the  helicopter  fuselage,  could  at  certain  frequencies  in  the  HF  band 
be  the  dominant  mode,  depending  on  the  mounting  and  orientation  of  the  loop  antenna. 
A  short  study,  restricted  to  the  2-15  MHz  band  which  is  used  for  HF  communications  over 
distances  less  than  about  600  km,  was  therefore  undertaken  to  test  this  hypothesis.  It  was 
not  the  purpose  of  the  study  to  validate  NEC2  for  wire-grid  models  yet  again.  Consequently 
no  attempt  was  made  to  ’tune’  the  numerical  model. 

2.  DESCRIPTION  OF  THE  NUMERICAL  AND  SCALE  MODELS 

In  order  to  investigate  the  radiation  patterns  of  a  loop  antenna  on  a  helicopter,  a  wire  grid 
structure,  which  can  best  be  described  as  a  generic  model  of  a  helicopter,  was  implemented 
for  use  in  the  Numerical  Electromagnetics  Code  NEC2  [6].  The  numerical  model  was 
based  on  top  and  side  view  scale  model  diagrams  of  an  SA330B  Puma  helicopter  available 
in  a  l/32nd  scale  kit  (no.  PK-507)  produced  under  the  registered  trademark  ’Matchbox’  by 
Lesney  Products,  London,  England.  The  full  scale  dimensions  of  the  helicopter  are  given 
in  Table  1.  Figure  1  shows  line  drawings  of  the  wiregrid  model  and  a  sideview  of  the  Puma 
helicopter.  Also  shown  is  a  simple  loop/dipole  model  suggested  by  Burke  [7]  for 
frequencies  less  than  the  fundamental  electrical  resonance  frequency. 

All  wires  used  in  numerically  modelling  the  full  scale  helicopter  fuselage  were  kept  shorter 
than  1.5  m  with  a  diameter  of  0.2  m,  except  where  noted.  For  the  analysis  the  longer 
dimension  of  the  helicopter  was  chosen  to  lie  along  the  x-axis.  The  loop  was  modelled 
using  three  straight  sections,  each  with  a  diameter  of  0.03  m.  The  fourth  side  of  the  loop 
was  formed  by  part  of  the  tailboom,  which  was  modelled  as  a  rectangular  grid  structure. 
The  overall  dimension  of  the  loop  was  0.3  x  2.6  m  and  it  was  positioned  on  the  starboard 
side  below  the  tailboom.  The  longer  loop  dimension  was  parallel  to  the  tailboom.  The 
plane  of  the  loop  was  at  an  angle  of  about  30  degrees  to  the  x-z  or  <p  =  0  degrees  plane, 
as  defined  by  the  rectangular  or  spherical  coordinate  system  in  the  IEEE  standard  test 
procedures  for  antennas  [8].  This  offset  in  the  plane  of  the  loop  was  decided  on  after 
looking  at  possible  means  of  attaching  the  loop  to  the  airframe.  The  feedpoint  was  in  the 
centre  of  the  short  loop  section  indicated  in  Figure  1(b).  The  shortest  segment  of  length 
0.06  m  occurred  in  the  feed  region,  giving  a  worst  case  segment  length  to  diameter  ratio  of 
2.  Except  for  the  number  of  main  rotor  blades  this  generic  model  in  general  resembles  that 
of  the  Sea  King  used  by  Kubina  in  [3].  The  total  number  of  wires  and  segments  used  were 
206  and  291  respectively.  The  numerical  modelling  guidelines  of  [3]  and  [4]  were  used. 

Initially  a  single  precision  version  of  NEC2,  restricted  to  a  maximum  of  300  wire  segments, 
was  used  in  order  to  speed  up  computer  throughput.  This  version,  however,  predicted 
physically  unacceptable  negative  radiation  resistances  at  frequencies  of  2  to  4  MHz  as  well 
as  different  null  depths  for  conducting  and  perfectly  conducting  wire  grids.  Furthermore, 
scaling  of  the  geometry  of  the  full-scale  numerical  model  at  HF  to  dimensions  appropriate 
for  VHF  model  measurements  by  using  the  GS  data  card  in  NEC2,  did  not  give  the  same 


^In  a  recent  paper  [11]  Cox  and  Vongas  draw  attention  to  the  fact  that  on  a  helicopter  a  loop  antenna  shunt 
feeds  the  fuselage/electric  dipole,  and  that  the  helicopter/loop  combination  radiates  in  both  loop  and  dipole  modes 
for  frequencies  less  than  the  lowest  electrical  resonance  frequency  of  the  helicopter. 
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results  as  the  unsealed  version.  The  problems  were  referred  to  Burke  who  ascribed  them 
to  the  use  of  electrically  short  segments  at  the  source  [7].  At  this  stage  of  the  investigation, 
an  improved  operating  system,  which  was  reputed  to  use  more  accurate  algorithms  for 
evaluating  functions,  was  installed  on  the  computer  used.  Comparison  of  outputs  run  under 
the  two  different  operating  systems,  showed  that  there  were  considerable  differences  (6  to 
10  dB)  in  null  depths  and  only  marginal  differences  in  peak  values  for  various  predicted 
radiation  patterns.  A  double  precision  version  of  NEC2,  which  also  resolved  the  problem 
of  negative  radiation  resistances,  was  therefore  implemented. 

The  numerical  modelling  was  independently  checked  by  a  colleague  on  a  VAX  computer 
[9],  but  with  smaller  grid  spacing  and  shorter  wirelengths.  This  model  used  880  elements. 
The  predicted  gain  results  at  5  MHz  were,  for  practical  purposes,  identical  to  those  obtained 
using  the  single  precision  version  of  NEC2  for  the  model  considered  in  this  paper.  Two 
cases  were  considered  in  the  independent  study  [9].  These  were  with  a  pair  of  the  main 
rotor  blades  aligned  with  the  fuselage  and  then  rotated  to  make  an  angle  of  45  degrees  with 
the  ’axis’  of  the  fuselage.  The  differences  in  the  predictions  for  these  two  cases  were 
negligibly  small  at  5  MHz,  suggesting  that  rotor  modulation  of  the  radiation  pattern  could 
be  small  at  frequencies  below  the  fundamental  resonance  frequency.  Although  this 
independent  analysis  indicated  that  the  proposed  numerical  modelling  was  probably 
adequate  for  the  purposes  of  this  study,  it  was  decided  to  validate  the  predicted  results  using 
a  scale  model. 

A  l/22nd  scale  model  of  the  helicopter,  based  on  the  kit  data  referred  to  above,  was 
manufactured  from  wood.  The  wooden  fuselage  was  covered  with  strips  of  thin  brass 
shimstock.  All  overlapping  joints  were  soldered  to  provide  electrical  continuity  in  any 
direction  across  the  interface  between  strips.  This  does  not  necessarily  model  the  expected 
surface  conductivity  of  the  helicopter  particularly  well  because  of  the  use  of  riveted  panels 
and  composite  materials  and  the  fact  that  the  metal  conductivities  cannot  always  be  scaled. 
No  allowance  was  made  in  the  model  for  doors  or  other  apertures.  The  loop  antenna  on 
the  model  was  fed  via  a  coaxial  cable  connection  on  the  inside  of  the  model.  The  rear 
rotors  as  well  as  the  horizontal  rear  stabiliser  were  made  of  metal.  All  provided  good 
electrical  contact  with  the  fuselage.  One  of  the  problems  which  arose  at  this  stage  was 
whether  to  model  the  main  rotors  using  carbon  fibre  composite  or  metal.  Examination  of 
a  rotorblade  of  a  Puma  at  a  South  African  Air  Force  open  day  showed  that  the  leading 
edge  was  protected  by  a  metal  strip  which  was  electrically  bonded  to  the  rotor  boss  and  the 
fuselage.  It  was  concluded  that  it  would  be  appropriate  to  treat  the  main  rotors  as 
conductors  as  had  been  done  in  the  numerical  modelling.  This  is  essentially  the  same 
approach  as  used  by  Owen  [4]. 

The  numerical  modelling  frequencies  for  the  22:1  scale  model  of  the  helicopter  ranged  from 
44  to  330  MHz  in  22  MHz  steps.  As  a  check  on  the  accuracy  of  the  computations,  gain 
calculations  were  made  at  5  degree  increments  in  azimuth  and  elevation  and  the  average 
gain  computed.  For  perfectly  conducting  wires  the  average  gain  in  free  space  should  be  1.0 
[6].  The  number  of  segments  of  the  wires  in  the  feed  region  and  other  wires  connected  to 
these  were  varied  in  order  to  obtain  the  best  overall  average  gain  over  the  frequency  range. 
The  position  of  the  feed  segment  was  also  changed  to  ensure  the  best  overall  average  gain. 
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The  lowest  segment  length  to  diameter  ratio  of  about  2  was  in  the  feed  region.  The 
extended  thin-wire  kernel  option  of  NEC2  was  tried  but  gave  no  real  improvement  in  the 
average  gain  predicted  over  the  frequency  range.  In  general  the  guidelines  for  numerical 
modelling  given  by  [3]  and  [4]  were  followed.  The  segment  lengths  of  all  wires  connected 
to  the  shorter  loop  sections  were  made  approximately  the  same  length  as  the  segments  used 
in  the  loop  sections.  The  diameters  of  these  wires  were  made  0.08  m  in  order  to  limit  the 
change  in  diameter  between  connected  wires  to  a  factor  of  about  2.5.  In  this  case  this  was 
found  to  give  satisfactory  results. 

Figure  2  shows  the  predicted  average  gains,  normalisation  gain  factors  and  directivities, 
corrected  for  the  deviation  in  the  average  gain  from  the  expected  value  of  1.0,  for  the 
numerical  model  of  the  22:1  helicopter  scale  model  using  perfect  and  brass  conductors.  The 
range  of  modelling  frequencies  from  44  to  330  MHz  corresponds  to  a  full  scale  frequency 
range  of  2  to  15  MHz.  Application  of  the  average  gain  correction  for  the  brass  wire  model 
gives  a  corrected  directivity  which  agrees  to  within  0.01  dB  with  that  shown  for  a  perfect 
conductor  in  Figure  2(a).  Below  200  MHz  the  directivity  is  approximately  that  expected  for 
a  dipole.  It  was  concluded  that  the  numerical  model  would  be  adequate  to  test  the 
hypothesis  that  the  dominant  radiation  mode  for  frequencies  less  than  the  fundamental 
resonance  frequency,  is  the  longitudinal  dipole  mode. 

A  limited  set  of  radiation  pattern  measurements  was  made  on  the  22:1  scale  model  at 
frequencies  of  110  and  220  MHz,  corresponding  to  full-scale  frequencies  of  5  and  10  MHz 
respectively.  All  measurements  were  performed  at  the  National  Antenna  Test  Facility  at 
Paardefontein,  north  of  Pretoria,  South  Africa  [10].  The  choice  of  frequencies  was  largely 
dictated  by  the  level  of  interference  in  the  FM  and  TV  broadcast  bands  at  the  site.  The 
source  antenna  was  a  horizontally  polarised  ground  arrayed  log-periodic  array.  A  10  dB 
attenuator  was  inserted  at  the  feedpoint  of  the  loop  antenna  in  order  to  provide  a  nearly 
constant  load  for  the  receiver.  This  restricted  the  overall  sensitivity  of  the  system  and  may 
have  influenced  the  measured  results.  The  orientation  of  the  model  was  adjusted  depending 
on  the  cut  and  required  field  component  relative  to  the  model’s  coordinate  system. 

3.  COMPARISON  BETWEEN  PREDICTED  AND  MEASURED  RESULTS 

The  maximum  path  length  on  the  helicopter  which  could  support  an  electrical  dipole-like 
resonance  was  estimated  to  be  about  19.2  m.  This  is  measured  along  the  length  of  a  rotor 
blade,  down  the  drive  shaft,  back  along  the  tailboom  and  out  to  the  end  of  the  rear  control 
surface.  For  an  electrical  dipole  this  corresponds  to  a  fundamental  resonance  frequency  of 
7.8  MHz,  corresponding  to  a  scaled  frequency  of  172  MHz.  For  the  helicopter  this 
frequency  may  be  reduced  because  of  the  diameter  to  length  ratio  of  the  fuselage  itself 
compared  with  a  thin  dipole,  and  possible  capacitive  loading  effects  of  the  main  and  tail 
rotors.  The  selected  measurement  frequencies  of  110  and  220  MHz  straddle  the 
fundamental  resonance  frequency  of  the  model. 

Because  of  the  difference  in  the  implementations  of  the  feedpoints  in  the  numerical  and 
physical  scale  models,  it  is  expected  that  a  comparison  of  the  predicted  and  measured  input 
impedances  for  the  loop  antenna  will  at  best  be  of  a  qualitative  nature.  The  feed  for  the 


56 


numerical  model  was  in  the  centre  of  the  short  segment  marked  as  ’feedpoinf  in  Figure 
1(b).  In  the  case  of  the  scaled  physical  model  this  loop  segment  was  fed  at  the  junction 
point  at  the  fuselage.  Figure  3  compares  the  predicted  and  measured  input  impedances  for 
the  physical  scale  model.  Only  half  of  the  Smith  chart  is  shown.  The  measured  results 
indicate  that  the  real  part  of  the  input  impedance  increases  from  1  Ohm  to  about  60  Ohms 
over  a  range  of  44  to  242  MHz.  Over  the  same  range  the  inductive  part  increases  from 
about  4  to  340  Ohms.  Figure  4  compares  the  measured  input  impedances  for  the  loop 
antenna  with  and  without  the  main  and  tail  rotors.  The  effect  of  removing  the  rotors  is 
relatively  small,  suggesting  that  rotating  rotors  will  not  modulate  the  loop’s  input  impedance 
by  much.  One  can  see  what  appears  to  be  a  relatively  small  effect  in  the  measured  input 
impedance  for  the  model  due  to  the  fundamental  resonance  frequency  at  176  MHz.  In  the 
case  of  the  fuselage  with  no  rotors,  the  probable  effect  of  a  slightly  higher  fundamental 
resonance  frequency  of  220  to  242  MHz  can  be  seen.  The  length  of  the  full  scale  fuselage 
excluding  rotors  is  about  13.8  metres,  corresponding  to  a  scaled  (model)  halfwave  resonance 
frequency  of  239  MHz. 

Ideally  an  examination  of  the  relative  amplitudes  of  the  theta-  and  phi-  components  of  the 
electric  field  in  the  x-y  plane  will  show  the  difference  in  the  relative  importance  of  the  loop 
and  electrical  dipole  fields,  E(0)  and  E(0)  respectively  in  this  case.  In  retrospect  it  would 
have  been  advisable  to  mount  the  loop  antenna  symmetrically  with  respect  to  the  tailboom 
and  below  it  in  the  x-z  plane,  instead  of  at  an  angle  of  30  degrees,  in  order  to  enhance  the 
differences  in  the  fields. 

For  all  helicopter  radiation  patterns  shown  in  the  figures  one  pair  of  rotors  was  aligned  with 
the  length  of  the  fuselage.  Also  shown  in  Figures  5,  6,  8  and  9  are  the  expected  normalised 
radiation  patterns  for  a  dipole  placed  along  the  ’axis’  of  the  fuselage.  Because  of  time 
limitations  on  the  use  of  the  test  range  no  measurements  were  made  in  the  vertical  y-z 
plane  orthogonal  to  the  fuselage  ’axis’.  The  alignment  of  the  rotors  in  the  helicopter 
schematic  in  the  figures  are  for  illustrative  purposes  only.  The  use  of  symbols  to  indicate 
orientation  are  as  follows:  F = forward,  S = starboard,  P = port,  T = top,  B = bottom  and  A = aft. 

Figures  5  and  6  show  normalised  comparisons  of  measured  and  predicted  radiation  patterns 
for  a  model  frequency  of  110  MHz,  corresponding  to  a  full-scale  frequency  of  5  MHz,  for 
x-y  and  x-z  planes  respectively.  In  Figure  5  the  model’s  measured  and  predicted  results  in 
the  horizontal  x-y  plane  at  1 10  MHz  (5  MHz  full-scale)  for  the  E(0)  and  E(0)  fields  are  in 
reasonable  agreement.  Except  for  the  depth  of  the  nulls,  the  E(0)  fields  resemble  the 
dipole  radiation  pattern.  The  E(0)  fields  can  be  attributed  mainly  to  the  loop  antenna  with 
null  filling  at  broadside  being  caused  by  vertical  currents  on  the  fuselage.  Figure  6  shows 
the  model’s  measured  and  predicted  radiation  patterns  at  110  MHz  for  the  E(0)  and  E(0) 
fields  in  the  vertical  x-z  plane.  In  this  case  the  E(0)  field  can  reasonably  be  expected  to  be 
a  combination  of  the  loop  and  dipole  modes.  The  predicted  and  measured  E(0)  patterns 
are  again  in  reasonable  agreement.  However,  they  tend  to  resemble  only  the  pattern  of  an 
electrical  dipole  aligned  with  the  fuselage  ’axis’,  with  little  contribution  evident  for  the  loop 
antenna.  Comparison  of  the  measured  and  predicted  E(0)  fields  is  however  poor.  Some 
of  the  differences  for  this  field  component  are  probably  due  to  the  differences  in  the  feed 
arrangements  for  the  numerical  and  physical  models  referred  to  earlier.  The  lack  of 
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measuring  sensitivity  and  stray  fields  due  to  scattering  from  coaxial  cables  may  also  have 
influenced  the  measurement  of  the  E(0)  component.  The  predicted  E(0)  field  resembles 
that  of  a  small  loop  antenna  in  the  x-y  plane.  This  may  be  a  residual  effect  of  the  offset 
loop  antenna  which  is  mainly  in  a  vertical  plane. 

A  residual  electrical  dipole  field,  believed  to  be  attributable  to  the  feed  arrangement  for 
the  numerical  model,  can  be  seen  in  the  E(0)  field  illustrated  in  Figure  7  which  shows  the 
predicted  fields  in  the  y-z  plane.  The  nulls  are  approximately  aligned  with  the  axis  of  the 
short  feed  section.  The  E(<p)  component  is  nearly  omni-directional,  with  a  variation  of  less 
than  1  dB,  as  would  be  expected  for  a  predominantly  electrical  dipole-like  field. 

Figures  8  and  9  compare  the  predicted  and  measured  E(0)  and  E(0)  fields  for  the  model 
at  220  MHz  (10  Mhz  full  scale)  for  the  x-y  and  x-z  planes  respectively.  It  should  be  noted 
that  the  rotor  blades  on  the  full  scale  helicopter  have  a  diameter  of  15  m,  corresponding 
to  a  resonant  halfwave  dipole  of  10  MHz.  In  the  case  of  the  x-y  plane,  shown  in  Figure  8, 
the  agreement  is  relatively  good  for  the  E(0)  component.  The  measured  E(0)  component 
once  more  tends  to  resemble  the  dipole  field,  but  compares  relatively  poorly  with  the 
predicted  field.  Figure  9  illustrates  the  case  for  the  x-z  plane.  Except  for  the  depth  of  the 
mill  the  agreement  between  the  measured  and  predicted  E(0)  fields  is  reasonable.  It  is 
believed  that  the  main  rotor  blade  resonance  influenced  the  predicted  depth  and  position 
of  the  null.  The  effect  of  changing  the  frequency  to  215  MHz  is  also  shown  for  comparison. 
The  E(0)  components  resemble  each  other  superficially.  Differences  are  ascribed  to 
problems  encountered  during  the  measurements  due  to  low  sensitivity  and  scattering  from 
coaxial  cables.  The  position  and  depth  of  null  of  the  predicted  E(0)  component  in  the  x-y 
plane  as  shown  in  Figure  8  was  found  to  be  not  nearly  as  sensitive  to  changes  in  frequency 
compared  to  the  E(0)  field  shown  in  Figure  9. 

The  results  illustrated  in  Figures  5  and  6  were  typical  of  predictions  up  to  about  176  MHz 
which  was  near  the  fundamental  electrical  resonance  frequency  of  the  model  helicopter. 
However,  it  was  found  that  as  the  frequency  for  the  scale  model  was  decreased  to  44  MHz 
(2  MHz  foil  scale),  the  predicted  contribution  of  the  loop  antenna  to  the  overall  radiation 
pattern  increased  steadily  until  the  loop’s  contribution  to  the  radiated  field  compared  with 
that  of  the  dipole,  as  shown  in  Figure  10.  Also  shown  are  the  fields  computed  for  the 
simple  electrical  dipole/loop  antenna  combination  of  Figure  1(c)  as  suggested  by  Burke  [7]. 
The  dimensions  of  the  electrical  dipole/loop  combination  are  similar  to  the  overall 
helicopter  and  loop  antenna  combination.  The  wire  diameters  were  0.02m.  The  null-filling 
evident  in  the  predicted  radiation  pattern  for  the  dipole  mode  for  the  helicopter  model  is 
ascribed  to  the  asymmetric  placement  of  the  loop  antenna,  with  resultant  circumferential 
currents  on  the  fuselage  in  the  x-y  plane.  Above  44  MHz  for  both  the  helicopter  and  simple 
case,  the  radiated  field  due  to  the  electrical  dipole  mode  exceeded  that  due  to  the  loop. 
In  both  cases  at  110  MHz  (5  MHz  full  scale)  the  radiated  field  due  to  the  electrical  dipole 
mode  was  predicted  to  be  about  15  dB  greater  than  that  due  to  the  loop.  At  220  MHz  (10 
MHz  full  scale)  for  the  simple  model,  the  electric  field  due  to  the  loop  is  predicted  to  be 
more  than  25  dB  less  than  that  due  to  the  electrical  dipole. 

Examination  of  the  predicted  radiation  patterns  for  the  model  helicopter  between  the 
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fundamental  resonance  frequency  of  about  176  and  330  MHz  (8  MHz  and  15  MHz  full 
scale),  shows  that  the  electrical  dipole  field  component  becomes  more  distorted.  Examples 
of  this  effect  are  shown  in  Figures  11  to  13  for  the  principal  planes  for  the  model  at  330 
MHz  (15  MHz  full  scale).  Figure  12  shows  a  similar  residual  electrical  dipole  effect, 
possibly  due  to  the  feed  segment  and  also  observed  in  Figure  7.  It  is  believed  that  the 
distortion  of  the  electrical  dipole  field  is  due  to  the  increasing  relative  importance  of  various 
higher  order  current  modes  discussed  in  [5]  as  the  frequency  is  increased.  It  is,  however, 
clear  that  there  is  still  a  substantial  contribution  from  the  electrical  dipole  mode  of 
radiation  for  the  helicopter. 

4.  DISCUSSION 

On  the  basis  of  reasonable  agreement  between  measured  and  predicted  results  for  the 
largest  electrical  field  components,  it  is  concluded  that  the  dominant  radiation  mode  for  a 
helicopter-borne  loop  antenna  lying  in  a  plane  containing  the  longest  dimension  of  the 
fuselage,  tends  to  be  electrical  dipole-like  at  frequencies  below  the  fundamental  resonance 
frequency  but  above  some  frequency  at  the  low  end  of  the  HF  band.  For  the  full  scale 
example  considered  here  this  low  end  frequency  is  of  the  order  of  2  MHz.  A  strong 
electrical  dipole-like  component  persists  up  to  15  MHz  for  the  example  considered. 
However,  as  frequencies  increase  above  the  fundamental  electrical  resonance,  the  effect  of 
higher  order  current  modes  on  the  fuselage  and  main  rotors  on  the  radiation  pattern 
become  increasingly  important.  The  resultant  pattern  is  thus  a  composite  of  the  electrical 
dipole-like  field,  fields  due  to  the  higher  order  current  modes  and  any  residual  effect  due 
to  the  loop  antenna  itself. 

’Nap  of  the  earth’  (NOE)  HF  communications  is  restricted  to  the  lower  half  of  the  HF  band 
by  ionospheric  propagation  and  require  high  take-off  angles.  Most  helicopters  and  aircraft 
will  have  fundamental  electrical  resonances  at  frequencies  in  the  HF  band  because  of  their 
physical  size.  For  this  reason  it  is  important  that  aircraft  antenna  engineers  concerned  with 
NOE  systems  consider  the  entire  helicopter  as  the  radiating  structure.  It  is  suggested  that 
it  could  be  worthwhile  examining  ways  in  which  the  dipole  mode  could  be  directly 
stimulated  in  order  to  improve  efficiency  and  matching. 
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Table  1. 


Dimensions  of  the  Puma  SA330  helicopter  used  in  the  numerical  and  scale 
models. 


(1) 

diameter  of  main  rotor  (4  blades) 

15  m 

(2) 

overall  length  (rotors  turning) 

18.15  m 

(3) 

height 

5.14  m 

(4) 

fuselage  length  including  rear  rotors 

14.82  m 

(scaled  from  kit  data) 
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Feedpoint 


Figures  (a)  and  (b),  isometric  and  side  views  of  the  wiregrid  model  of  a 
SA330  PUMA  helicopter;  (c),  simple  low  HF  loop/dipole  model;  (d)  line 
drawing  of  the  helicopter. 
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Average  gains,  predicted  total  gain  normalisation  factors  and  corrected 
directivities  for  the  helicopter  model  with  (a)  perfectly  conducting  and  (b) 
brass  wires. 


[•] 


gure  3.  Comparison  of  predicted  and  measured  antenna  input  impedances  for  the 
complete  helicopter  scale  model.  Only  the  upper  half  of  the  Smith  chart  is 
shown.  (Normalisation  -  50  ohms) 
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Complete  model 


All  rotors  removed 
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Comparison  of  measured  antenna  input  impedances  for  the  complete 
helicopter  scale  model  and  the  fuselage  after  removal  of  all  rotors.  Only  the 
upper  half  of  the  Smith  chart  is  shown.  (Normalisation  -  50  ohms) 


Figure  5.  Comparison  of  normalised  measured  and  predicted  radiation  patterns  at  110 
MHz  in  the  x-y  plane.  Predicted  patterns  normalised  to  E(0)  =  1.1  dB. 
Normalised  dipole  pattern  shown  for  reference. 
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Figure  6. 


Comparison  of  normalised  measured  and  predicted  radiation  patterns  at  110 
MHz  in  the  x-z  plane.  Predicted  patterns  normalised  to  E(0)  =  1.3  dB. 
Normalised  dipole  pattern  shown  for  reference. 


Figure  7. 


Predicted  fields  for  the  helicopter  model  in  the  y-z  plane  at  110  MHz. 
Patterns  normalised  to  E(0)  =  1.5  dB. 
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Figure  8.  Comparison  of  normalised  measured  and  predicted  radiation  patterns  at  220 
MHz  in  the  x-y  plane.  Predicted  patterns  normalised  to  E(0)  =  0.6  dB. 
Normalised  dipole  pattern  shown  for  reference. 


Figure  9.  Comparison  of  normalised  measured  and  predicted  radiation  patterns  at  220 
MHz  in  the  x-z  plane.  Predicted  patterns  normalised  to  £(6)  =  1.1  dB. 
Normalised  dipole  pattern  shown  for  reference. 
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Figure  10.  Comparison  of  predicted  radiation  patterns  at  44  MHz  for  the  scaled 
helicopter  borne  loop  antenna  and  the  simple  loop/ dipole  model,  using  brass 
conductors  for  both  cases.  Helicopter  patterns  normalised  to  E(0)  =  -16.5 
dB.  Loop/dipole  patterns  normalised  to  E(0)  =  -18.0  dB.  The  loop  is  in  the 
x-z  plane.  Patterns  are  for  the  x-y  plane. 


Figure  11.  Predicted  gains  for  the  helicopter  model  in  the  x-z  plane  at  330  MHz.  Brass 
conductors. 
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Figure  12.  Predicted  gains  for  the  helicopter  model  in  the  y-z  plane  at  330  MHz.  Brass 
conductors. 
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Figure  13.  Predicted  gains  for  the  helicopter  model  in  the  x-y  plane  at  330  MHz.  Brass 
conductors. 
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Numerical  analysis  of  reffector  antennas  uses  a  discrete  approximation  of  the  radi¬ 
ation  integral.  The  calculation  replaces  the  actual  reflector  surface  with  a  triangular 
facet  representation.  The  physical  optics  current  is  then  approximated  within  each 
facet.  This  paper  provides  analytical  details  of  the  method  based  on  the  assump¬ 
tion  of  a  constant  magnitude  and  linear-phase  approximation  of  the  physical  optics 
current.  Example  calculations  are  provided  for  parabolic,  elliptical,  and  shaped 
subreflectors.  The  computed  results  are  compared  with  calculations  made  using  a 
constant-phase  approximation.  These  results  show  that  the  linear-phase  approxi¬ 
mation  is  a  significant  improvement  over  the  constant-phase  approximation  in  that 
the  solution  converges  over  a  larger  angular  region  of  space.  This  improvement  can 
significantly  reduce  storage  requirements  and  possibly  execution  speed. 


I.  Introduction 

One  of  the  simplest  numerical  techniques  for  electromagnetic  scattering  analysis  is  based  on  a  discrete 
approximation  of  the  physical  optics  (PO)  radiation  integral.  The  general  modeling  technique  is  similar  to 
that  employed  by  Rao  et  al.  [1]  for  the  moment  method  solution  of  electromagnetic  scattering  by  surfaces 
of  arbitrary  shape.  In  this  paper,  we  apply  the  methodology  to  the  problem  of  shaped  reflector  antenna 
analysis.  This  calculation  comprises  two  distinct  approximations.  First,  the  actual  reflector  surface  is 
replaced  by  a  triangular  facet  representation  so  that  the  reflector  resembles  a  geodesic  dome.  One  then 
makes  an  analytic  approximation  of  the  PO  current  within  each  of  the  facets.  Upon  evaluating  the  PO 
integral  locally  over  each  facet,  the  radiation  integral  reduces  to  a  summation  over  the  collection  facets  that 
represent  the  surface. 

Several  years  ago,  a  computer  program  was  developed  at  the  Jet  Propulsion  Laboratory  (JPL)  utilizing 
the  assumption  of  constant  magnitude  and  phase  of  the  PO  current  within  each  facet.  This  program  has 
proven  to  be  surprisingly  robust  and  useful  for  the  analysis  of  relatively  small  shaped  reflectors,  particularly 
when  the  near  fleld  is  desired  and  surface  derivatives  are  not  known.  It  is  natural  to  inquire  whether  a  more 
sophisticated  approximation  of  the  PO  surface  current  will  yield  more  accurate  results  or  permit  the  use 
of  larger  facets.  In  this  paper,  a  linear-phase  approximation  of  the  surface  currents  is  made.  Within  each 
triangular  region,  the  resulting  integral  is  written  as  the  two-dimensional  Fourier  transform  of  the  projected 
triangle.  This  triangular  shape  function  can  be  integrated  in  closed  form  [2]  and  the  complete  PO  integral  is 
then  a  summation  of  these  transforms.  Significantly,  other  authors  have  developed  more  general  techniques 
for  performing  the  required  integration  [3,4,5],  which  could  be  very  useful  for  future  refinements. 

In  what  follows,  the  explicit  details  of  the  analysis  are  provided  along  with  example  calculations  of 
scattering  from  parabolic,  elliptical,  and  shaped  reflector  surfaces.  For  a  given  size  of  triangular  facet,  two 
general  trends  emerge  from  the  calculations.  First,  the  linear-phase  approximation  takes  about  three  times 
longer  to  compute  a  field  point  than  does  the  constant-phase  approximation.  Second,  there  is  an  angular 
region  of  space  over  which  the  solution  is  valid,  and  this  angular  region  is  significantly  larger  with  the  linear- 
phase  approximation  than  with  the  constant-phase  approximation.  Clearly,  a  trade-off  situation  exists  here. 
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Since  the  surface  geometry  and  PO  currents  must  be  stored  in  memory,  the  linear-phase  approximation 
offers  an  advantage  in  terms  of  storage  because  fewer  triangular  facets  are  needed  to  reach  convergence  over 
a  specified  angular  region.  On  the  other  hand,  to  claim  a  speed  advantage,  this  method  must  reduce  the 
required  number  of  facets  by  at  least  one-third.  This  reduction  will  only  be  possible  for  reflectors  in  which 
relatively  large  regions  can  be  adequately  approximated  by  a  uniformly  illuminated  planar  surface.  It  has 
been  found  that  this  is  the  case  for  relatively  large  reflector  antennas. 


11.  Analytical  Details 

The  PO  radiation  integral  over  the  reflector  surface  E  can  be  expressed  [6] 

1  /*  /  1  \  ^  ^-jhR 


(1) 


in  which  r  designates  the  field  point,  r'  the  source  point,  R  =  |r  — r'|  is  the  distance  between  them,  and  R  — 
(r  —  r')/i?  is  a  unit  vector.  The  PO  surface  current  on  the  subreflector  surface  is  expressed 

J,(r')  =  2iixH,(r')  (2) 


For  the  purpose  of  analysis,  the  true  surface  E  is  replaced  by  a  contiguous  set  of  TV-plane  triangular  facets. 
These  facets,  denoted  A*,  are  chosen  to  be  roughly  equal  in  size  with  their  vertices  on  the  surface  E.  Figure 
1  shows  a  typical  facet  and  its  projection  onto  the  x-y  plane.  Let  represent  the  centroid  of  each 

triangle  where  the  subscript  i  =  1,  •  —  ,  TV  is  associated  with  a  triangle.  Then,  the  field  obtained  by  replacing 
the  true  surface  E  by  the  triangular  facet  approximation  is 


H(r) 
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R 
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(3) 


In  Eq.  (3),  J  is  now  the  equivalent  surface  current  evaluated  on  the  triangular  facets.  Since  the  triangles 
are  small,  it  is  expected  that  R  and  R  do  not  vary  appreciably  over  the  area  of  a  given  facet.  Thus,  let  R^ 
and  Ri  be  the  value  obtained  at  the  centroid  {xi^yi^zi)  of  each  facet  and  approximate  Eq.  (3)  by 


H(r)  =  - 


Air 


Ki  X  Ti(r) 


(4) 


(5) 


Assume  that  the  necessary  transformations  have  been  performed  so  that  the  incident  field  is  given  in 
terms  of  the  reflector  coordinate  system.  Then 

Ji(r')  =2n,- X  H,(r')  (6) 
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Next,  assume  that  the  incident  field  can  be  represented  by  a  function  of  the  form 


H,  =h,(r'i) 


e-jkr. 

4Trr,i 


(7) 


where  is  the  distance  to  the  source  point.  Also,  r^-  and  denote  the  vectors  r'  and  evaluated  at  the 
centroid  of  the  ith.  triangular  facet  and  are  regarded  as  constants.  Then,  Eq.  (5)  can  be  written 

Ti(r)  =  (8) 

27r  RiVg{  JAi 

To  simplify  the  form  of  the  integration,  the  surface  Jacobian  is  introduced  within  each  triangular  facet 
For  a  planar  surface  Zi  =  fi{x,y),  a  normal  is  given  by 


where 


and  a  unit  normal  is  given  by 


Nj  =  -xfri  -  yfyi  +  z 


f  =dfi  ,  ^dfi 


Ni 
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This  permits  the  explicit  evaluation  of  the  Jacobian  as 


Ja,  =1  N, 


1/2 


fxi  +  fyi  +  1 


(11) 


Making  use  of  the  Jacobian  then  allows  Eq.  (8)  to  be  rewritten  as 


Ti{v) 
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2'KRiTgi 


Ja, 


(12) 


in  which  A'-  represents  the  area  of  the  zth  triangular  facet  projected  onto  the  z  ~  plane.  Now,  make  a 
Taylor-series  expansion  of  the  exponent  in  Eq.  (12).  Retaining  only  the  first-order  terms,  one  can  formally 
write 


R{x,  y)  +  rs{x,  y)  =  -{ai  -  UiX  -  viy) 


(13) 


in  which  a^,  Uj,  and  Vi  are  constants.  This  approximation  corresponds  to  a  far-field  approximation  on  the 
fth  triangle.  With  this  approximation,  Eq.  (12)  reduces  to 
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It  may  now  be  observed  that  this  integral  is  the  two-dimensional  Fourier  transform  of  the  ith  projected 
triangle  A'-,  expressed  as 


S{u,v)=  I  (15) 

In  order  to  explicitly  evaluate  the  constants  in  Eq.  (13),  note  that  the  equation  of  a  plane  can  be  expressed 
as 


Z  —  {x  4"  (y  Vi^fyi  4“ 


This  can  be  used  to  obtain 
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Placing  the  result  of  Eq.  (16)  into  Eq.  (14),  and  recalling  Eqs.  (6)  and  (7),  yields 
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This  is  the  final  form  of  the  linear-phase  approximation  over  each  triangular  facet.  This  expression  can  be 
used  in  Eq.  (4)  to  compute  the  radiation  integral  once  the  Fourier  transform  of  a  triangular  shape  function 
S{u^v)  is  known.  Fortunately,  this  transform  can  be  computed  in  closed  form  [2]  from  the  expression 
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(20) 


in  which  {xn,yn)  are  the  coordinates  of  the  triangle  vertices  numbered  in  a  clockwise  direction.  The  slope 
of  the  nth  side  (between  corners  n  and  n  -1-  1)  is  given  by 


Vn+l  ~  Vn 
^n+1  Xn 


(21) 


Some  attention  must  be  given  to  the  following  special  cases.  First,  if  u  0,  the  transform  reduces  to 

the  formula  for  the  area  of  a  triangle 


5(0,0)  =  -i 


xiiy2  -  2/3)  +  X2{y3  -  yi)  +  -  2/2) 


(22) 
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Next,  \i  u/v 


-Pn,  then 
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III.  Numerical  Results 

A  FORTRAN  subroutine  was  written  to  perform  the  linear-phase  calculations  indicated  above.  Test  cases 
were  run  for  parabolas,  ellipses,  and  shaped  subreflectors,  and  the  results  were  compared  with  calculations 
that  use  a  constant  magnitude  and  phase  approximation  on  the  triangular  facets.  A  focused  parabola  is 
neither  an  interesting  nor  a  challenging  case  for  the  algorithm,  since  the  phase  variation  over  the  facet  is 
small.  As  a  simple  test  case,  the  far-field  pattern  and  gain  of  a  1,  OOOA-diameter  parabolic  reflector  with  a 
focal  length  of  F  =  400A  was  calculated.  The  reflector  is  illuminated  by  a  linearly  polarized  horn  with  a  cos  0 
pattern  function.  Figure  2  compares  the  linear-  and  constant-phase  approximation  for  a  roughly  equally 
spaced  80-by-80  rectangular  grid  of  points  divided  into  triangles  over  the  reflector  surface  (approximately 
10,000  triangles).  The  running  time  on  a  Cray  X-MP/18  was  less  than  one  minute.  Convergence  was 
checked  by  increasing  the  number  of  triangles  until  the  computed  solution  did  not  change  appreciably  over 
the  angular  region  of  interest.  It  has  been  previously  demonstrated  [7,8,9]  that,  once  sufficient  triangles 
to  converge  the  solution  have  been  utilized,  the  results  of  the  constant-phase  algorithm  are  valid,  so  only 
comparisons  of  the  two  techniques  are  presented. 

A  more  interesting  example  is  the  ellipse  shown  in  Fig.  3.  The  projected  aperture  of  the  ellipse  is  about 
3  m,  illumination  function  is  a  cos^^  ^  pattern  function  (22.3-dB  gain),  and  the  frequency  is  31.4  GHz. 
The  ellipse  is  about  350 A  along  the  major  axis.  Figure  4  compares  the  constant-phase  approximation  for 
different  grid  densities  of  approximately  4,000,  10,000,  and  23,000  triangles,  and  illustrates  a  general  trend 
of  the  method,  i.e.,  depending  on  the  size  of  the  triangles,  there  is  an  angular  limit  over  which  the  solution 
is  valid.  Figure  5  compares  the  linear-phase  approximation  with  the  constant-phase  approximation  for  the 
4,000-triangle  case  and  demonstrates  that  the  angular  range  is  larger  with  the  linear-phase  approximation. 

A  third  example  is  the  shaped  subreflector  shown  in  Fig.  6.  The  diameter  is  3.42  m  (135  in.),  and  it  is 
fed  with  a  cos^^^  ^  pattern  function  (29.7-dB  gain).  Figure  7  compares  the  results  of  a  4,000-  and  10,000- 
triangle  grid  constant-phase  approximation  with  a  4,000-triangle  linear-phase  approximation.  The  frequency 
of  operation  is  2.3  GHz,  hence,  the  subreflector  is  about  26A  in  diameter.  The  10,000-triangle  constant  phase 
is  the  converged  result  and  the  4, 000- triangle  linear  case  gives  the  same  result.  A  very  good  approximation 
is  also  obtained  with  a  1,400-triangle  grid  for  the  linear  case,  but  no  meaningful  results  are  obtained  with 
the  constant-phase  case.  Figure  8  gives  the  linear-phase  result  for  31.4  GHz  (360A  subreflector)  using  23,000 
triangles.  No  meaningful  result  is  obtained  for  the  equivalent  constant-phase  case. 

One  final  example  is  given  by  a  beam- waveguide  system  that  has  recently  been  built  at  JPL.  The  measure¬ 
ment  setup  consists  of  a  22-dB-gain  feedhorn  that  is  used  to  illuminate  a  be  am- waveguide  system  consisting 
of  a  pair  of  parabolic  reflectors.  The  mirrors  are  arranged  to  replicate  the  input  feed  pattern  at  the  focal 
point  of  the  second  dish.  Details  of  the  geometry  are  given  elsewhere  [10].  A  calculation  of  this  system 
was  made  using  the  triangular  facet  PO  technique  by  first  computing  the  near-field  scattering  from  the  first 
reflector  and  using  these  values  to  compute  the  PO  current  on  the  second  reflector.  Subsequently,  the  tri¬ 
angular  facet  program  is  used  a  second  time  to  compute  the  field  radiated  by  the  second  parabolic  reflector. 
In  Fig.  9,  the  results  of  this  calculation  are  compared  with  experimental  measured  data  taken  at  X-band. 
As  can  be  seen,  the  computed  results  compare  favorably  with  the  measured  data. 

Most  of  the  examples  given  are  for  large  reflectors  to  illustrate  the  robust  character  of  the  technique. 
For  smaller  reflectors  (<  100 A),  meaningful  results  can  be  obtained  on  a  typical  desktop  computer  in  a 
reasonable  time. 
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IV.  Conclusions 

The  triangular  facet  approximation  technique  provides  a  simple  and  flexible  method  of  analysis  for  re¬ 
flector  antennas.  We  have  found  that  the  linear-phase  approximation  is  valid  over  a  larger  angular  region 
than  is  the  constant-phase  approximation.  In  applying  this  method  to  fairly  large  (100  to  lOOOA)  reflectors, 
the  linear-phase  approximation  provides  a  significant  reduction  in  the  computer  storage  requirements  and 
can  reduce  computation  time  in  some  cases. 
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Fig.  1.  Reflector  analysis  coordinate  systems  and  a  typical 
triangular  facet. 
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Fig.  4.  Ellipse  example:  constant-phase  approximation  for  offset  plane. 
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Fig.  5.  Ellipse  example:  constant  versus  linear  phase  for  offset  plane. 


SHAPED 

SUBREFLECTOR 


Fig.  6.  Shaped  subreflector  geometry. 
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Fig.  7.  Shaped  subreflector  example  for  H-plane  at  2.3  GHz. 
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Abstract 

Wire-grid  modelling  of  continuous  surfaces  for  structures  in  three  dimensions  is  a  tedious  and  time 
consuming  process.  This  paper  describes  a  simplified  automatic  mesh  generator  that  converts  a  large 
class  of  three-dimensional  structures  into  appropriate  wire-grid  models.  The  output  of  the  generator 
can  readily  be  used  as  the  input  to  wire  antenna  Moment  Method  codes.  It  is  however  designed 
specifically  for  use  with  the  Numerical  Electromagnetics  Code  (NEC). 


1.  Introduction 

Users  of  Moment  Method  codes  often  require  the  generation  of  wire-grid 
models  for  three-dimensional  structures  such  as  ships,  vehicles,  aircraft,  etc. 
Producing  these  models  manually  usually  takes  a  lot  of  time  and  effort.  Even 
when  working  from  the  manufacturers  CAD  drawings  the  user  often  needs  to 
eliminate  unnecessary  details  and  parts  of  the  structure.  Although  some 
sophisticated  computer  software  has  been  developed  for  this  purpose  [1],  it  is  not 
generally  available. 

The  Mesh  Generator  developed  in  this  work  can  produce  wire-grid  models 
for  these  3-D  structures  according  to  the  required  specifications.  The  structure 
needs  to  be  represented  by  flat  surfaces  of  quadrilateral  or  triangular  shape.  The 
grids  produced  by  the  generator  are  generally  quasi-orthogonal  wires  except  near 
one  of  the  vertices  in  triangular  regions;  as  will  be  shown  later.  Limited  modelling 
errors  are  checked  for  by  the  generator;  although  not  as  extensively  as  described 
in  [2]. 

The  output  of  the  generator  can  be  used  as  a  NEC  [3]  input  data  file.  External 
graphics  facilities  need  to  be  used  to  view  the  generated  wire-grid. 


2.  Mesh  Generation  Guidelines 

Wire-grid  modelling  techniques  have  been  studied  extensively  by  many 
workers.  It  is  not  the  intention  of  this  paper  to  give  any  recommendations  in  this 
regard.  However,  because  of  their  relevance,  mesh  generation  guidelines  will  be 
briefly  reviewed  here.  When  modelling  a  solid  surface  with  a  wire-grid,  three  main 
factors  need  to  be  considered: 
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•  Type  of  the  wire-grid 

•  Cell  size  of  the  grid 

•  Radii  of  grid  wires 

The  most  widely  used  type  of  wire-grid  is  the  rectangular  one  because  it  was 
found  to  be  satisfactory  in  most  cases  and  it  is  easier  to  generate  than  other  types. 
Also,  the  rectangular  grid  makes  it  easier  to  decide  on  the  values  of  wire  radii 
using  simple  formulas.  A  triangular  mesh  on  the  other  hand  has  the  advantage 
of  simulating  more  realistic  currents  [4,  p.117].  Other  types  of  meshes  have  also 
been  used.  The  mesh  generator  described  in  this  paper  can  produce  rectangular 
and  partly  triangular  grids. 

The  grid  cell  size  can  either  be  specified  in  terms  of  its  area,  as  in  [4],  or  in 
terms  of  the  spacing  between  adjacent  elements.  Because  the  rectangular  mesh 
is  the  one  that  is  most  commonly  used,  it  is  easier  to  use  the  second  method.  The 
element  spacing  (size)  will  depend  on  the  type  of  structure  being  modelled  and 
on  the  parameters  of  interest.  It  must  be  as  small  as  possible,  up  to  the  limitations 
of  the  MM  code  used,  if  these  are  near  field  parameters  such  as  input  impedance. 
The  grid  element  size  can  range  from  X/5,  as  suggested  in  [5],  to  as  little  as  V1500 
as  used  in  [6]  near  the  feed  point  of  the  antenna.  The  generally  accepted  value  for 
element  size  is  X/10. 

As  for  the  wire  radii  of  the  grid,  two  basic  methods  have  been  discussed  in 
the  literature.  The  first  specifies  the  wire  radius  of  the  model  as  a  fraction  of  the 
wavelength.  Studies  that  have  used  this  method,  such  as  in  [7],  set  the  wire  radius 
empirically  at  0.005X  but  do  not  give  any  reasoning  for  this  choice.  The  second 
method  is  much  more  widely  used  and  relates  the  overall  surface  area  of  the 
wires.  As  ,  to  the  area  of  the  solid  surface  that  is  being  modelled,  Ar .  As  has  been 
described  in  [8],  modellers  have  recommended  several  criteria  for  this  relation. 
The  following  table  summarises  those  most  commonly  used:* 


Modeller 

Criterion 

Moore  &  Pizer 

As  >  Ar  up  to  5Ar 

Lee,  et  al 

Ludwig 

Elliott  &  Mcbride 

A^  —  Ar  of  orw  side  of  a  closed  surface 

As  —  Ar  of  both  sides  of  an  open  surface 

Burke  &  Poggio 

As  of  wires  in  the  direction  of  each  current 
component  =  Ar .  So,  for  a  rectangular 
mesh,  As  =  2Ar 

'  Complete  references  for  these  criteria  are  found  in  [8]. 
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It  is  evident  that  most  of  these  criteria  recommend  the  "twice  surface  area 
rule",  i.e.  As  =  2Ar  .  If  we  define  the  ratio  between  As  and  Ar  to  be  the  overall 
area  factor,  AF  ,  then  it  follows  that  AF  =  2  in  the  general  case.  It  is  however 
useful  to  keep  AF  as  a  variable  set  by  the  user  according  to  what  is  required.  This 
mesh  generator  uses  the  given  area  factor  to  calculate  the  wire  radius  for  the  mesh 
of  the  modelled  surface. 


3.  Description 

The  program  for  this  mesh  generator  is  written  in  standard  FORTRAN  77. 
It  can  run  on  a  mainframe  or  a  microcomputer.  The  generator  is  used  to  model 
structures  in  the  following  manner: 

1.  The  structure  is  divided  into  quadrilateral  or  triangular  flat  regions  of  any 
orientation  in  space. 

2.  Knowing  the  corners  co-ordinates  of  each  individual  region  and  the  number 
of  sections  required  on  its  sides,  the  region  is  converted  into  a  wire-grid 
accordingly. 

3.  The  wire  radius  for  each  region  is  calculated  according  to  the  specified  area 
factor  as  described  below. 

4.  Similar  elements  at  the  junctions  of  regions  sides  are  then  eliminated. 

5.  Parts  of  the  structure  that  do  not  form  flat  surfaces  must  be  added  manually. 

A  quadrilateral  region  is  defined  by  the  co-ordinates  of  its  four  corners  Cl 
to  C4.  There  is  no  restriction  on  the  shape  of  the  quadrilateral  region  or  its 
orientation  in  space  as  long  as  it  is  flat.  In  the  special  case  of  the  triangular 
region,  the  co-ordinates  of  two  consecutive  corners  are  the  same.  The  region  is 
divided  into  grid  elements  by  specifying  the  number  of  elements  on  two  adjacent 
sides  (N1  and  N2  in  Fig.  1). 

The  other  two  opposite  sides  are  assumed  to  have  the  same  number  of 
elements  as  those  opposite  regardless  of  their  lengths.  For  triangular  regions, 
there  are  only  three  corners,  so  one  of  the  opposite  sides  will  simply  be  a  vertex. 
To  establish  a  reference,  this  type  of  region  is  defined  as  one  in  which  corners  3 
and  4  have  the  same  co-ordinates.  Consequently,  there  are  three  alternatives  for 
the  mesh  of  a  triangular  region  as  shown  in  the  three  schematics  of  Fig.  2.  By 
changing  the  numbering  order  of  the  region's  three  corners,  the  required  mesh 
can  be  generated.  N1  and  N2  were  kept  the  same  at  4  and  3  respectively,  in  all 
the  three  combinations  in  order  to  illustrate  the  point.  In  practice  they  can  be 
changed  to  suit  the  model. 

The  generator  will  then  produce  all  elements  that  make  up  the  wire-grid  of 
the  region.  The  sequence  of  ordering  the  elements  is  shown  in  Fig.  1  for  the  first 
four  elements  (LI  to  L4). 
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The  next  step  is  to  calculate  the  radius  of  wires  in  the  region  using  the  given 
area  factor.  The  total  surface  area  of  the  wires  As  is  first  calculated  by  the 
following  formula: 

k 

/=1 

where 

//  :  The  length  of  the  ith  element 

a  :  The  radius  of  the  region's  wires 
k  :  The  number  of  elements  in  the  region 


The  region's  area  Ar  is  then  calculated  by  dividing  the  region  into  two 
triangles  whose  areas  are  found  using  Heron's  formula  [9,  p.7]  as  follows: 

Af.  =  Ai-\-  A2 
^1  = 


^12  +  ^23  +  ^31 


^2  =  —  ■^34)(‘^2  “  ■^4l)(‘^2  ~  -^b) 


-^34  +  -^41  +  -^13 
2 


In  the  above  expressions,  Lf^n  refers  to  the  length  of  the  region's  side 
connecting  corners  m  and  n.  The  generator  will  then  determine  the  wire  radius 
a,  which  is  the  same  for  all  the  elements  of  the  region,  using  the  expression: 

AFxAr 

k 

/=1 


The  above  process  is  repeated  for  each  region  in  the  structure.  Since  common 
elements  exist  at  the  junctions  of  regions'  sides,  all  are  scanned  and  the  duplicated 
elements  are  eliminated.  The  wire  radius  of  the  elements  at  a  junction  is 
determined  using  the  parameters  of  the  region  listed  first  in  the  input  file. 

At  this  stage,  most  of  the  structure  will  have  been  converted  automatically 
to  a  wire-grid  of  the  required  specification.  The  user  only  has  to  add  those  parts 
of  the  structure  that  could  not  be  represented  by  flat  surfaces.  These  for  example, 
may  include  curved  surfaces  or  single  wires. 

It  is  important  to  note  that  connected  regions  should  have  the  same  number 
of  elements  for  the  sides  at  the  junctions.  It  is  also  advisable  that  they  have 
similar  mesh  densities  so  that  any  change  in  stepped  wire  radius  between  two 
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connected  regions  is  made  as  gradual  as  possible.  Different  mesh  densities  for 
different  parts  of  the  structure  can  be  achieved  by  making  a  gradual  change  in 
the  mesh  densities  of  connected  regions.  After  all,  the  limitations  of  the  Moment 
Method  code  that  is  used  will  determine  the  stepped  radius  factor  that  it  can 
handle. 

In  addition,  the  generator  produces  some  useful  information  for  the  different 
regions  of  the  structure  and  checks  for  some  of  the  common  modelling  errors.  The 
following  parameters  are  produced  for  each  region  in  the  structure: 

•  Segment  length  4  in  wavelengths 

•  Wire  radius  a  in  wavelengths 

•  Ratio  of  segment  length  to  wire  radius  4/a 

The  user  can  specify  the  error  threshold  for  these  parameters  or  use  the  default 
values  which  are  set  according  to  the  suggested  guidelines  in  [2].  These  are: 

4  >  /i/5 
XI a  <  30 
Is!  a  <  2 

The  threshold  for  the  4/«  ratio  is  set  at  2  assuming  that  the  extended  thin-wire 
kernel  in  NEC  will  not  be  used. 

Note  that  the  information  given  for  individual  regions  are  for  all  the  elements 
in  the  region  including  those  that  will  be  eliminated  when  the  regions  are  joined 
together. 


4.  Inputs  and  Outputs  of  the  Generator 

The  input  data  of  the  generator  is  read  in  free  format  from  a  file  structured 
as  follows: 


NREG  AF  FREQ 
XXI  XX2  XX3  XX4 
YYl  YY2  YY3  YY4 
ZZl  ZZ2  ZZ3  ZZ4 
N1  N2 

where 


FREQ 

NREG 

AF 

XXI  ...  XX4 
YYl  ...  YY4 
ZZl  ...  ZZ4 


Frequency  of  operation  in  Mega  Hertz. 

The  number  of  regions  on  the  structure. 

The  area  factor  for  all  the  regions. 

X  co-ordinates  of  the  four  corners  of  a  region  in  metres. 
Y  co-ordinates  of  the  four  corners  of  a  region  in  metres. 
Z  co-ordinates  of  the  four  corners  of  a  region  in  metres. 
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N1 

N2 


:  Number  of  elements  on  the  side  between  corners  1  and  2 
;  Number  of  elements  on  the  side  between  corners  2  and  3 


Data  for  all  the  regions  are  input  sequentially  in  the  same  manner  one  after 
the  other  without  leaving  spaces.  Note  that  the  area  factor,  AF  ,  is  specified  only 
once  at  the  beginning  of  the  file  for  the  entire  structure. 

Each  element  generated  is  simply  a  one-segment  wire.  Two  output  files  are 
produced  by  the  generator: 

Geometry  output  file:  which  is  formatted  as  a  standard  NEC  geometry  file  that 
can  be  processed  using  the  NEEDS  package  [10]  or  a  mainframe  version  of  NEC. 
The  file  format  is  accepted  by  both  IGUANA  [10]  and  NECPLOT  [3]  graphics 
routines. 

Information  file:  which  includes  the  following: 

•  The  area  factor,  number  of  regions,  and  total  number  of  elements  for  the 
structure. 

•  Ar,  As,  a,  and  the  total  length  of  elements  for  each  region  given  in  metres 
squared  and  metres. 

•  The  modelling  parameters  mentioned  in  the  previous  section,  grouped  for 
each  region  separately,  and  the  modelling  errors  if  found. 

It  was  suggested  by  an  anonymous  reviewer  that  one  option  to  this  generator 
would  be  to  specify  the  number  of  segments  per  wavelength  between  corners  of 
a  region  and  request  the  wire-grid  to  be  generated  accordingly.  While  this  option 
will  be  quite  useful  for  the  user,  its  implementation  would  require  substantial 
review  of  the  software. 


5.  Examples 

To  illustrate  the  usefulness  of  this  mesh  generator,  three  examples  will  be 
given.  The  first  is  a  rectangular  box  consisting  of  6  regions,  one  for  each  face,  and 
104  elements.  A  monopole  is  mounted  in  the  middle  of  its  top  side.  The  area 
factor  specified  for  the  box  is  2.0  and  the  frequency  of  operation  is  30  MHz.  The 
input  file  is  listed  below  for  illustration: 

6  2.0  30.00 

0.  2.  2.  0. 

0.  0.  0.  0. 

0.  0.  3.  3. 

2  3 

2.  2.  2.  2. 

0.  4.  4.  0. 

0.  0.  3.  3. 

4  3 
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0.  2.  2.  0. 
4.  4.  4.  4. 
0.  0.  3.  3. 

2  3 

0.  0.  0.  0. 
0.  4.  4.  0. 
0.  0.  3.  3. 
4  3 

0.  2.  2.  0. 
0.  0.  4.  4. 
0.  0.  0.  0. 
2  4 

0.  2.  2.  0. 


Fig.  (3)  Wire-grid  of  a  rectangular  box  consisting  of 
104  elements  and  a  monopole  on  top  of  it. 


Only  the  monopole  is  added  manually  after  the  wire-grid  has  been  generated. 
Fig.  3  is  an  isometric  view  of  this  simple  structure.  The  second  example  in  Fig. 
4,  is  a  pyramid  with  a  trapezoidal  base  where  triangular  regions  have  been  used. 
It  consists  of  5  regions  and  1 57  elements. 

The  last  example  is  more  practical  in  that  it  is  a  wire-grid  of  a  motor  vehicle. 
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Fig.  (4)  Wire-grid  of  a  pyramid  with  a  trapezoidal  base 
consisting  of  5  regions  and  157  elements. 


All  the  body  can  be  represented  by  flat  surfaces  except  the  pillars  that  support 
the  roof.  Because  the  number  of  elements  where  the  regions'  sides  join  needs  to 
be  the  same,  14  regions  were  used  to  represent  the  vehicle  body: 

1  Region  for  the  roof 
1  Region  for  the  bonnet  (hood) 

1  Region  for  the  boot  (trunk)  top 
1  Region  for  the  front  side 
1  Region  for  the  back  side 
3  Regions  for  the  left  side 
3  regions  for  the  right  side 
3  Regions  for  the  bottom. 

Fig.  5  shows  the  wire-grid  of  270  elements  that  was  produced  by  the  mesh 
generator.  In  order  to  get  the  final  model  of  the  vehicle,  the  following  steps  were 
done  manually: 


•  On  each  side  of  the  vehicle,  remove  the  two  vertical  elements  marked  x  in  Fig. 
5;  extend  the  wire  marked  y  to  the  left  so  that  it  meets  with  the  wire  junction 
J  in  the  adjacent  region;  and  join  the  top  of  these  two  regions  by  adding  a 
wire  between  points  pi  and  p2. 
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•  Add  the  wires  for  the  roof  pillars, 

•  Increase  the  number  of  segments  of  the  two  roof  elements  marked  z  from  1 
to  2  so  that  pillars  join  the  roof  at  segment  ends,  as  required  by  NEC. 

The  final  model  is  shown  in  Fig.  6. 


Fig.  (5)  Wire-grid  of  a  motor  vehicle  as  produced  by  the  mesh 
generator.  It  consists  of  14  regions  and  270  elements. 


Typical  execution  times  for  these  three  examples  are  shown  in  the  table  below. 
They  were  run  using  a  VAX  11/780  and  an  80386  IBM-PC  with  a  coprocessor 
using  the  FTN77  compiler.  Times  shown  are  those  taken  to  generate  the  NEC 
geometry  file  only  without  the  information  file. 


Example 

Run  Time  (sec.) 

VAX 

IBM-PC 

Box 

2.5 

Pyramid 

4.3 

Vehicle 

10.0 
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Fig.  (6)  Wire-grid  of  the  motor  vehicle  after  editing  the 
generated  model  and  adding  the  roof  pillars. 


6,  Availability  of  the  Software 

The  software  for  this  mesh  generator,  including  the  source  code  and  the  input 
and  output  files  for  the  examples  mentioned  above,  are  available  from  the  author. 
A  nominal  fee  of  £5.00  (or  an  international  money  order  for  US$10.00)  is 
requested  to  cover  duplication  and  postage.  The  software  will  be  provided  on  a 
3.5"  MS-DOS  720K  floppy  disk.  Hardcopy  supplementary  information  will  also 
be  provided. 
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THE  SHORT  FAT  DIPOLE 
Development  in  APL  of  a  MoM  Solution 

Douglas  B.  Miron 
EE  Dept.,  Box  2220 
South  Dakota  State  University 
Brookings,  SD  57007 

Abstract 

This  paper  describes  an  effort  to  find  the  best  profile  for 
a  figure-of-revolution,  center-fed,  electrically-small  dipole. 

It  includes  presentations  on  equation  development,  singularity 
treatment,  code  development,  verification,  and  performance. 

I.  Prologue 

The  work  described  in  this  paper  was  undertaken  to  find  a 
minimum-Q  shape  for  an  electrically-small  dipole  to  operate  over 
the  HF  band,  3-30  MHz.  The  starting  point  was  the  1969  paper  by 
Mautz  and  Harrington[l]  on  bodies  of  revolution.  Their  formula¬ 
tion  for  the  integral  equation  and  their  expansion  and  testing 
functions  were  used  without  scaling.  The  dipole  is  driven  only 
in  the  tangential  direction,  so  it  was  assumed  there  are  no 
0-components  to  the  current  density.  This  simplified  the 
problem  considerably.  There  are  numerous  ways  to  handle  the 
potential— function  singularity.  After  trying  a  couple  of  others, 
subtraction  of  the  electrostatic  potential  from  the  radiation 
potential  function[2]  was  chosen.  The  electrostatic  potential 
function  is  integratable  in  0  to  a  known  special  function,  and 
that  function  is  numerically  integratable  so  that  the  singularity 
is  expressed  in  the  impedance-matrix  elements  in  a  stable  way. 

This  paper  is  divided  into  the  following  sections: 

I.  Prologue 

II.  Geometry,  the  Basis  Functions,  and  Counting 

III.  The  Integral  Equation 

IV.  Conversion  to  Discrete  Form 

V.  Verification 

VI.  Impedance  and  Shape 

VII.  Development  Time  vs.  Execution  Time 
Appendix  APL  Syntax,  Symbols  and  Functions 

One  purpose  of  this  paper  is  to  demonstrate  that  programming 
in  an  interactive  array-processing  language  (APL)  makes  the  most 
efficient  use  of  professional  time.  The  Appendix  briefly  de¬ 
scribes  its  basic  syntax  and  lists  the  definitions  of  symbols  and 
functions  as  used  in  this  paper.  Specific  features  are  de¬ 
scribed  as  they  are  used  in  Sections  II  through  IV  .  In  Section 
V  the  work  is  tried  by  comparison  with  two  classic  results.  In 
Section  VI  a  search  for  a  shape  that  gives  an  antenna  Q  close  to 
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the  theoretical  minimum,  l/()8a)^  [3],  is  outlined.  In  Section 
VII  the  hardware-software  performance  issues  are  examined. 

II.  Geometry,  the  Basis  Functions,  and  Counting 

The  dipole  is  symmetric  and  center-fed.  It  must  fit  in  a 
2  m  cube  and  be  easy  to  manufacture.  It  would  be  a  good  thing 
then,  if  its  profile  were  made  of  connected  straight-line  seg¬ 
ments.  Figure  1  shows  an  example  of  a  simple  shape,  along  with 
some  of  the  coordinates  to  be  used  in  the  calculations.  All 
prospective  shapes  are  to  be  figures  of  revolution  about  the  z 
axis,  with  the  feed  gap  at  z=0.  An  x  axis  is  necessary  from 
which  to  define  the  rotational  angle  0.  Cylindrical  radius  p 
is  used  both  for  the  shape  description  and  in  the  calculations. 
The  fundamental  coordinate  for  integration  is  the  tangential 
distance  t,  which  runs  from  t=0  at  z=-l  m  to  t=t^  at  z=+l  m, 
following  the  antenna's  profile. 


Figure  1.  A  representative  Figure  2.  x  and  z  projections 

dipole  as  a  figure  of  revo-  of  a  tangential  vector.  Xp  is 

lution  about  the  z  axis.  parallel  to  the  x  axis. 

One  must  know  the  projection  of  the  current  density  vector 
at  one  place  onto  the  surface  at  another  place,  to  calculate  the 
vector  magnetic  potential  and  surface  E.  Since  the  applied 
field  is  in  the  z  direction,  it  is  assumed  there  are  no  0 
components  of  current  or  electric  intensity.  Figure  2  shows  the 
components  of  projected  onto  the  z  and  x  axes. 
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(1) 


Jp  =  J^sinO',  =  Jocose',  =  JjSine'cos^' 

For  the  purpose  of  doing  a  dot  product  with  a  unit  tangential 
vector  at  the  observation  point,  the  last  projection  should  be  by 
way  of  cos (0-0'),  but  the  rotational  symmetry  allows  0=0. 

The  triangle  function  was  used  as  the  basic  element  of  the 
current  expansion  and  weighting  function  series.  The  triangle, 
which  is  shown  in  Figure  3,  was  divided  by  p(t)  to  make  the 
complete  scalar  function.  Since  the  integrals  were  discretized 
it  was  easier  to  denote  the  triangle  by  a  subscript  and  deal  with 
its  relation  to  t  through  their  counters,  which  is  displayed  in 
Figure  4.  However,  for  the  sake  of  the  functional  notation  used 
in  the  next  section,  one  may  say  that 


Figure  3.  (a)  Triangle  function  and  a  piecewise-constant 

approximation.  (b)  The  derivative  of  a  triangle 
function. 


Tj  (t) 

W,.(t)  =  UtWi(t)  = 

where  u^  is  a  unit  vector,  W,-  are  the  vector  basis  functions, 
and  their  divergences  are 


i  ^ 


The  dipole  can  be  described  by  assigning  vectors  of  p  and  z 
values  at  the  corners  in  the  profile  of  the  half-dipole.  For 
example 


RHDt-0.01  10  0  ZD«-0  0.9  1 
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/v  A  /v  / 

/  \  /  'v  /  /  \  / 

/  \  /  \  /  \  /  ^  ' 


X  X  A  A 


/  /  \  /  \  /  \  /  \ 
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K=  1  23456789  10 


Figure  4.  Relationship  between  the  counter  for  t, 
tk=(k-J)Dt,  and  the  triangle  counter. 

Tj  spans  from  k=2i-l  to  k=2i+2. 

describes  a  dipole  which  has  a  radius  of  0.01  m  at  the  feedpoint, 
cones  out  to  a  1  in  radius  at  z=0 . 9  m  and  cones  in  to  a  point  at 
z=l  m.  The  full  dipole  description  is  formed  by  reversing  ZD  and 
RHD,  negating  the  reversed  ZD  values,  and  attaching  them  to  the 
originals. 

ZV«-(-(l)ZD)  ,0,ZD  (4) 

RHV<-(<t)RHD)  ,RHD[1]  ,RHD  (5) 

The  extra  point  in  the  middle  is  a  convenience  so  the  first  z 
point  doesn't  have  to  be  zero,  and  a  straight  pipe  may  be  insert¬ 
ed  before  the  shape  takes  off. 

Figure  5  shows  the  basis  for  a  scheme  to  generate  the 
coordinate  variables  as  functions  of  t.  Since  the  profile  is  a 
straight  line  between  corners,  the  values  of  p  and  z  on  a 
segment  are  displacements  from  their  corner  values  which  are 
sin0  and  cos0  times  the  t  displacement  from  TC[K] ,  where  K  is  a 
corner  counter.  The  nub  of  the  scheme  is  to  generate  the  TC 
vector,  and  then  write  a  function  to  find  the  index  of  the  last 
corner  for  a  given  t.  The  increase  in  t  between  corners  is  the 
length  of  the  profile  segment. 

TC[K+l]«-TC[K]  +  (  (  (RHV[K+1]-RHV[K] )  *2)  +  (ZV[K+1] -ZV[K]  )  *2)  *0. 5 

This  expression  could  be  put  in  a  loop  and  stepped  through  to 
produce  the  vector  TC.  However,  APL  has  a  cumulative  sum  opera¬ 
tor,  and  there  are  a  couple  of  simple  ways  to  generate 

the  pairwise— difference  of  a  vector.  Here  is  one  as  a  defined 
function. 
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Figure  5.  Illustration  of  a  piecewise-linear  antenna  profile 

section,  z  and  t  increase  to  the  right,  p  increases 
upward,  a  and  0  are  positive  ccw. 


7DIFF[D]V 
[0]  Z«-DIFF  X 

[1]  Z*-(liX)--liX  (6) 

Now  the  expression  for  TC  is  written 

TC*-0,+\(  (  (DIFF  RHV)  *2)  +  (DIFF  ZV)*2)*0.5  (7) 

For  a  single  value  of  t,  one  way  to  find  the  index  of  TC  just 
left  of  t  is  to  compare  the  value  with  TC.  TC<T  will  produce  a 
Boolean  vector  with  1  for  each  corner  up  to  the  nearest  one  on 
the  left  of  T,  and  0  for  the  rest.  The  index  value  for  the 
corner  wanted  is  the  sum  of  these  Is,  +/TC<T.  To  do  this  for  a 
vector  of  t  values,  the  outer  product  operator  can  be  used. 
TCo.<,T  takes  each  value  of  TC  and  compares  it  to  all  the 
values  in  T,  producing  a  row  in  the  result.  The  result  has  as 
many  columns  as  there  are  elements  in  T  and  as  many  rows  as  there 
are  elements  in  TC.  The  ' next  to  T  makes  a  vector  of  a 
single-valued  T.  To  get  the  index  values  now  it  is  necessary  to 
sum  down  the  columns  of  the  Boolean  matrix.  These  steps  are  put 
in  a  defined  function  called  REGION. 
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VREGIONEDIV 
[0]  S<-REGION  T 

[1]  S<-+/TC®.<,T  (8) 

In  this  function  TC  is  used  as  a  global  variable,  T  and  S  are 
local  variables  through  which  data  is  passed  into  and  out  of  the 
function. 

For  the  purpose  of  numerical  integration,  the  profile  of  the 
antenna,  which  is  the  range  of  t,  is  divided  into  increments  of 
A  and  Dt  as  in  Figure  3 .  The  center  of  the  dipole  has  to  be  a 
center  point  for  an  expansion  function  and  a  boundary  for  a  A 
segment.  Let  NP  be  the  number  of  A  segments  in  the  half-dipole. 

DEL*-TC[pTC]^2xNP  and  DT*-DELt2  (9) 

NP  is  also  the  number  of  unique  current  coefficients,  so  it  is 
quite  important  later  on. 

Finally,  ^2it/K  (BT<-o2t LAMBDA)  and  ^  (BS^-BTxBT)  are 
important  values  in  the  calculations  because  /3  measures  the 
intervals  in  wavelength.  LAMBDA  serves  as  a  global  variable  to 
hold  K. 

III.  The  Integral  Equation 
IIIA  Field  Theory 

The  antenna  is  a  closed,  perfectly-conducting  (initial 
assumption)  surface  so  that  the  currents  flowing  on  it  and  the 
applied  voltage  must  produce  tangential  E  fields  that  add  to 
zero.  The  field  due  to  the  surface  currents  is 

=  -Jo'"Jo1jw/^(t')g(R)-7g(R)^^^^^^)p'd0  (10) 

I 

The  distance  function  is 


R  =  -2pp ' cos (0-0' )  +  ( z-z ' )  (11) 

and  the  propagation  function  is 
-ifiR 

g(R)  =  W 
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The  5-gap  source  model  was  used  so  the  incident  E  is 


E’(r)  =  u,V„5(z)  (13) 

and  the  boundary  condition  requires  that  (13)  equal  the  negative 
of  the  tangential  component  of  (10) .  This  would  be  an  ill-condi¬ 
tioned  integral  equation  on  several  counts,  but  MoM  smooths  it 
out  somewhat  by  multiplying  both  sides  by  a  weighting,  or  averag¬ 
ing,  function  and  then  integrating  the  products  over  the  surface 
of  the  antenna.  Using  a  set  of  vector  weighting  functions  and 
the  vector  dot  product  (inner  product)  also  picks  out  the  tangen¬ 
tial  component  on  each  side.  Representing  a  weighting  function 
by 

W(t)  =  UtW(t)  (14) 

the  smoothed  equation  becomes 


J^"'J^’'w.E’pd0dt  = 


h(t).j(tMg(R) 


+W(t).7g(R) 


V'-j(t') 


p'd0'dt'pd0dt 


Again,  from  symmetry  there  are  no  0-components  to  the  sources  or 
fields,  so  the  left  side  integrates  to 


27rp(V2)W(V2)V„ 

The  first  dot  product  on  the  right  can  be  represented  by  the  unit 
vectors  involved, 

Uj'Uj,  =  cos0cos9'+  sinOsinG'cos  (0-0' )  (16) 

from  Figure  2  and  the  discussion  around  equation  (1).  The  second 
dot  product  involves  the  gradient  of  g(R),  which  is  already  a 
singular  function.  Although  a  second-order  singularity  can  be 
handled  [4,5],  Harrington  and  Mautz  [1]  used  the  fact  that  the 
body  is  a  closed  surface  to  replace  the  integral  of  W(t)'Vg(R) 
in  an  integration  by  parts  with  -g(R)V”W  (t). 

R  and  the  unit-vector  dot  product  involve  cos (0-0'). 

These  are  the  only  rotational  terms,  so  they  were  collected  and 
the  angle  integrals  done  first.  Because  the  functions  really 
depend  on  the  difference  between  the  two  angles,  if  one  inte¬ 
grates  with  respect  to  one  angle,  the  dependence  on  the  other 
disappears.  The  0-dependent  terms  can  be  collected  in  two 
integral  functions  as 
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(17) 


“  j)3R 

°o  =  2irJ^'g(R)d,>  = 

and  likewise 

-  -jiSR 

Gi  =  J^COS0^^  d0  (18) 

Collecting  all  these  substitutions  brought  the  integral  equation 
to 


p(^)W(^)Vo 

j60;3 


=  J*^'"J^"'{W(t)  J(t' )  (cos0(t)cos9(t')Go+ 


sin0(t)sine(tMGi)  -  GqV-W  (t)  V'-J(t' )//g^}PP'dt'dt  (19) 


IIIB  The  Integrals  in  0 

Gq  and  Gi  hold  the  integrations  in  0,  which  must  be  done 
before  the  outer  integrals  in  t  and  t'.  They  also  hold  the  1/R 
singularity.  There  must  be  some  advantage  to  the  rotational 
symmetry  besides  the  simplicity  of  the  setup,  perhaps  it  would  be 
in  the  evaluation  of  these  integrals.  At  first,  writing 

jSR=^(R-R,)+iSR,,  rI=P^+(z-z')^ 

was  tried,  R^  being  the  distance  to  the  axis  under  the  source 
point.  For  /3(R-Ro)<l  a  few  terms  in  the  exponential  series 
were  used,  which  then  led  to  terms  involving  elliptic  integrals 
of  the  three  kinds.  While  this  was  educational,  and  some  errors 
in  [6]  were  found,  the  result  was  rather  complex  and  slow.  An 
interesting  discovery  was  that,  away  from  the  singularity,  it  is 
faster  to  compute  the  integral  by  the  piecewise-constant  series 
approximation  than  to  use  the  polynomial-and-log  approximation 
from  [7].  The  speed  crossover  was  at  20  terms  in  the  series  for 
six-figure  accuracy. 

The  next  approach  tried  was  an  old  one  used  recently  by 
Simpson,  et.  al.  [2].  They  split  the  exponential  instead  of  /SR. 


T  e  1 
R  R  " 


(20) 


This  method  removes  the  singularity  from  the  propagation  function 
to  a  static  function  which  integrates  to  elliptic  integrals. 

These  two  parts  can  be  thought  of  as  the  circuit-element  part  and 
the  radiation  part,  concepts  which  are  explored  in  [2].  Define 
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~  Jq 


(21) 


^1E  “  Jr 


COS0  . 


-  J 


^1u.  “  J 


cos(^d0 


One  can  find  the  elliptic-function  version  of  Gq^  either  by 
a  succession  of  substitutions  or  from  [6]  as 


with  a  =  p^+p'^+  (z-z '  )^  f  b  —  2pp'  and  r‘^  — 

K(r)  is  the  complete  elliptic  integral  of  the  first  kind  whose 
defining  equation  is 


2  _  2b 


-7/2 

K(r)  =  J  -p=P 
Jl-r^s 


Using  the  same  substitutions  and  the  redundancy [8 ] 

sin^z  =  p  “  p(l-r^sin^z) 
provides 


E(r)  is  the  complete  elliptic  integral  of  the  second  kind, 
defined  by 


,ir/<!  I - - 

E(r)  =  4l^ 


sin^xdx 


No  closed-form  or  special-function  versions  for  Gq^ 
and  G^^  were  found,  but  the  integrands  are  smooth  and 
well-behaved  so  they  were  simply  approximated  by  the  same 
piecewise-constant  interval  integration  as  used  by  MoM. 
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me  The  Singularity 


The  development  in  IIIB  moved  from  the  static  potential 
function,  which  has  a  1/R  singularity  through  the  integration  in 
<p  to  the  elliptic  integrals.  [7]  gives  a  quite  accurate  approx¬ 
imation  which  explicitly  includes  a  logarithm  which  goes  to 
infinity  when  r-»l,  which  occurs  when  t'=t.  For  t-t'  small 
compared  to  p  the  argument  of  the  approximation  is  proportional 
to  t-t',  so  the  approximation  can  be  integrated  analytically. 
However,  there  are  other  functions  which  multiply  the  Gs,  and 
t-t'  is  not  necessarily  small  compared  to  p  over  the  interval  in 
the  discrete  approximation  to  the  integrals. 

The  basic  approach  by  MoM  to  solving  equation  (20)  is  to 
break  up  the  integrals  into  sums,  each  term  of  which  approximates 
the  integral  over  a  small  interval  by  the  product  of  the  interval 
length  and  the  integrand's  value  at  the  center  of  the  interval. 
This  is  the  piecewise-constant  approximation.  The  question  then 
is,  what  value  should  be  used  to  represent  the  G  functions  for 
intervals  in  which  t'-»t?  Working  the  charged-cylinder  electro¬ 
static  problem[9,  10]  showed  that  displacing  the  source  and 
observation  points  by  a  small  amount  compared  to  the  interval 
size  gave  a  result  independent  of  this  displacement,  once  it  was 
small  enough,  as  far  as  plots  of  the  charge  density  were  con¬ 
cerned.  In  the  present  case,  the  results  were  poor  for  small  NP, 
even  having  sign  errors  in  the  terminal  impedance.  Apparently 
the  singularity  was  over-represented.  The  next  best  choice 
seemed  to  be  to  use  the  average  value  of  each  G  function,  which 
means  numerically  integrating  them  over  the  interval  since  no 
analytical  integration  was  available.  Reference  [7]  provided 
some  formulas  and  values  for  functions  with  a  logarithmic  singu¬ 
larity.  The  numerical  approximation  has  the  form 

J^f  (x)  log(x)dx  =  i:w,f(x,)  (30) 

where  w,-  and  Xj  are  given  in  tables  for  various  numbers  of 
samples  per  interval.  The  form  needed  was 

IG  =  J*^G(y)dy  =  2j^G(y)dy  (31) 

since  G  is  even  in  y.  If  y=5x,  then  the  upper  limit  becomes  1, 


IG  =  2  J’g(<Sx)  5dx 

To  use  (30),  f(x)  =  G ( (Sx) /log (x)  , 


IG 


2(5Sw,- 


G(«Sx,) 

log(Xi) 


(32) 


(33) 
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The  problem  at  hand  was  to  integrate  over  a  Dt  interval  with 
respect  to  t'  when  t  is  at  the  interval  center,  y  corresponds  to 
t'-t,  6  to  Dt/2.  The  variables  in  the  G  functions  were  expressed 
as  functions  of  p,  p',  z  and  z',  so  they  were  re-expressed  in 
terms  of  t'-t,  p  and  sine(t)  .  This  was  made  possible  by  assuming 
t'  and  t  are  on  the  same  profile  segment.  The  geometry  of  this 
assumption  gives 

(p-p')^+(z-z^)^  =  (t'-t)^  (34) 

pf  =  p+(t'-t)sin9  (35) 

so  that 

a  =  (t'-t)^+2p(p+(t'-t)sine)  (36) 

b  =  2P (P+ (t'-t) sin9)  (37) 

a+b  =  (t'-t)^+4P(P+(t'-t)sin9)  (38) 

a-b  =  (t'-t)^  (39) 


™  ^  ^  a+b  a+b  '  ' 

m  is  the  most  convenient  variab.1  <=  in  which  to  express  the  ellip¬ 
tic  integral  approximations. 

IV  Conversion  to  Discrete  Form 

IVA  Using  MoM 

On  examining  equation  (20) ,  one  can  see  that  each  weighting 
function  is  multiplied  by  a  P  and  each  occurrence  of  the  current 
density  is  multiplied  by  a  P'.  Furthermore,  the  first  deriva¬ 
tive  of  both  functions  is  required.  It  makes  good  sense,  in 
hindsight,  that  Mautz  and  Harrington[ 1]  chose  to  represent  both 
functions  as  a  series  of  basis  functions  having  a  first  deriva¬ 
tive  and  being  a  simple  function  divided  by  P.  The  weighting 
functions  are 

W,.(t)  =  1-i-N  (41) 

The  surface  current  density  is 


_  a-b  _  (t^-t) 


j(t')  = 


One  can  show,  by  taking  the  textbook  approach  of  sketching  a 
differential-sized  box  on  the  antenna  surface  and  applying  the 
definition  of  divergence  as  a  Gauss's  Law  limit,  that 
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(43) 


l3(pW,-) 

P  3t 


V'-J  =  S  Jtk^W+rr 
k=i  ‘•'P'dt' 

There  can  easily  be  a  notation  and  counting  confusion  because  the 
value  of  t  at  the  middle  of  a  triangle  is  never  actually  used  in 
the  integral  computations.  It  is  simpler  to  express  the  inte¬ 
grals,  though,  using  this  value,  so  let  it  be  t,-  for  triangle  T,- . 
Then,  since  the  triangle  functions  are  nonzero  over  a  2A  domain 
the  integrals  are  shortened  up.  The  application  of  each 
weighting  function  to  equation  (20)  generates  a  separate 
equation,  the  full  sequence  produces  a  set  of  N  equations  for  the 
N  Jj  amplitudes. 


T,  (-^)V<, 
160^0 


S  Tj  (t)T,.(t')  (cosecose'Go+sinesine'G.) 

1C=1 


Go  dT:  (t)  dT|,(t')  1 

4  -4^4-  dt'dt,  i<i<N 


The  left-hand  sides  of  this  equation  sequence  are  zero  except 
for  i=NP  because  only  that  triangle  function  is  nonzero  at  t^/2 . 
The  results  of  the  double  integration  on  the  right-hand  sides  are 
functions  only  of  i  and  k,  so  they  can  be  symbolized  by  the 
standard  Zj,^. 


V,-  =  S  JtkZfk.  (46) 

k=1 

IVB  Preparing  the  Way;  Geometry 

The  function  values  needed  in  each  sine,  T,-,  etc.,  are 

used  over  and  over  again,  so  it  doesn't  pay  to  recompute  them 
each  time.  The  basic  strategy  was  to  precompute  them  as  vectors 
and  matrices,  and  call  the  needed  data  by  their  indices  during 
the  Z  element  calculation. 


APL  provides  a  function  to  generate  a  vector  of  index 
values,  Any  set  of  values  which  can  be  written  as  a 

function  of  an  index  can  be  generated  as  a  vector,  making  a  loop, 
with  its  repetitive  interpretation,  unnecessary.  One  can  gener¬ 
ate  a  function  of  any  number  of  indices  by  using  the  outer 
product  operator  with  the  individual  index  vectors  which  will 
produce  an  array  with  as  many  axes  as  there  are  index  variables. 
While  this  is  quick  to  execute,  large  temporary  outer  products 
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can  eat  up  storage  to  the  point  where  there  isn't  enough,  which 
will  generate  a  WS  FULL  (WS  for  workspace)  error  message. 

The  vectors  and  static  G  functions  are  calculated  in  a 
function  SGENERATE. 

[ 0 ]  SGENERATE ; K ; T ; THT ; L ; T1 ; N1 ; N2 

All  the  variables  in  the  header  are  local. 

[1]  Tl<-nTS  O  Nl<-2xNP  O  N2«-4xNP 

[2]  ZT«-ZA  T«-DTx'0.5+iN2  O  RHT«-RHO  T  O  THT«-TH  T  O  C0ST<-20THT 
O  SINT^-IOTHT 

Remember  that  the  interpreter  works  from  right  to  left  in  a 
statement.  The  statements  in  a  line  are  executed  starting  with 
the  leftmost  one.  The  first  statement  in  line  [2]  starts  by 
generating  a  vector  of  integers  from  1  to  N2,  then  adds  '0.5  to 
these  integers,  then  multiplies  them  by  DT  and  assigns  the  result 
to  T,  which  is  the  vector  of  t  values  used  in  all  subsequent 
functions.  The  T  is  also  an  input  argument  to  a  function  ZA 
which  calculates  the  corresponding  vector  of  z(t)  values,  ZT. 

The  next  two  statements  in  line  [2]  calculate  P(t)  and  0(t) 
vectors,  RHT  and  THT.  All  three  of  these  functions  use  the 
simple  segment  geometry  and  the  function  REGION  (6) .  The  next 
two  statements  in  line  [2]  produce  cos9  and  sinO  using  the 
circle  function  with  a  left  argument  to  indicate  which  of  the  15 
transcendental  functions  is  wanted. 

[3]  TRT«-(1  3  3  1)t4  O  DTRT<-(1  1  '1  ■l)^DEL 

These  are  vectors  to  hold  the  sample  values  of  the  triangle 
function  and  its  derivative. 

[4]  G0*-Gl^-(Nl,N2)p0  O  K<-1 

GO  and  G1  are  matrices  which  hold  the  static  potential  integrals 
G()£  and  G^g.  These  matrices  must  be  symmetric  about  their  main 
diagonal,  regardless  of  the  geometry  of  a  problem,  and  they  are 
also  symmetric  about  a  line  between  rows  N1  and  Nl+1  because  of 
the  geometric  symmetry  of  the  dipole.  Use  was  not  made  of  the 
diagonal  symmetry  because  of  the  complexity  of  the  indexing 
required.  The  row  symmetry  was  used. 

[5]  (iNl)  REALGREEN  lN2  O  k^l 

REALGREEN  is  a  function  which  calculates  the  values  for  the 
static  G  functions  and  puts  them  in  the  arrays  GO  and  Gl.  Its 
input  arguments  are  sets  of  row  and  column  indices.  The  workings 
of  REALGREEN  and  IG  will  be  discussed  below. 
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[6]  GL:]>-IG  RHT[K]  ,SINT[K]  O  GO  [K;K] «-L[  1]  O  G1[K;K]*-L[2] 

O  -»(N1>K*-K+1)/GL 

This  is  a  one-line  loop  to  replace  the  diagonal  elements  of  GO 
and  G1  with  the  average  values  of  Gq^  and  G^^  on  the  intervals 
where  they  are  singular.  The  calculation  on  these  intervals  in 
REALGREEN  were  prevented  from  overflow  by  replacing  zeros  of  m 
with  1E‘50.  It  would  have  been  more  efficient  to  send  the 
results  directly  to  GO  and  G1  inside  IG,  but  it  was  left  in  its 
development  form  which  gives  an  explicit  result. 

[7]  GO<-GO,  [l]0GO[Nl+l-i3;]  O  G1*-G1,  [  1] 0G1  [Nl+l-i3  ;  ] 

Not  all  the  possible  values  of  GO  and  G1  are  ever  used.  Again 
because  of  geometric  symmetry,  there  are  only  NP  unique  values  of 
Jji^,  so  only  NP  equations  are  needed.  The  expression 
Nl+1-13  generates  a  vector  Nl,  Nl-l,Nl-2  which  picks  these  three 
rows  out  of  the  matrix.  No  column  index  is  specified  so  the 
whole  row  is  taken.  This  three-row  block  is  then  reversed 
end-for-end,  and  attached  to  the  original  matrix  under  its 
columns,  making  the  result  a  2NP+3  by  4NP  array. 

[8]  (+/0  0  0  3600  60  1  1E‘3xdTS-T1)  'SEC,  STATIC  GO  AND  Gl.' 

This  line  calculates  and  displays  the  elapsed  time  in  seconds. 

The  mixture  of  numeric  and  character  data  in  a  vector  is  a 
feature  of  APL2 .  In  standard  APL  the  number  would  have  to  be 
converted  to  characters  to  be  displayed  in  the  same  vector. 

[9]  DSOUND  SCALE, [1.5] 100 

This  line  plays  an  ascending  C  major  scale  to  get  the  user's 
attention  when  the  function  is  finished  executing. 

IVC  Preparing  the  Way;  Potential  Integrals 


REALGREEN  is  so  named  because  its  results  are  Green's 
functions  of  a  sort. 

[0]  OX  REALGREEN  SX;A; B;SQAB;M;K;E 

OX  and  SX  stand  for  observation-point  indices  and  source-point 
indices,  respectively.  These  are  the  input  arguments,  the  other 
variables  listed  in  the  header  are  local  to  the  function  and 
named  to  correspond  to  variables  in  (25)  and  (28)  .  SQAB='Ja+b. 

[1]  B«-2xRHT[0X]<^  .xRHT[SX]  O  AB<-RHT [OX] . +RHT [SX] 

O  A«-ZT[OX]o.-ZT[SX] 
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In  each  of  the  statements  in  line  [1]  the  generalized  outer- 
product  operator  forms  a  matrix  of  all  combinations  of  the  left 
argument  elements  with  the  right  argument  elements,  with  the 
arithmetic  function  placed  between  them.  B  has  as  many  rows  as 
there  are  values  in  OX  and  as  many  columns  as  values  in  SX,  and 
each  element  of  B  is  2p  (t, )  p  (t,^)  =2pp ' .  Each  element  of  AB  is  p+p', 
and  each  element  of  A  is  z-z'. 

[2]  AB«- (AxA) +ABXAB  O  A<-OpO  O  50+ZERO  l-2xB^AB 

The  first  statement  forms  a  matrix  of  a+b  values.  The  second 
statement  reduces  A  to  an  empty  vector  to  reduce  storage.  The 
third  statement  calculates  an  array  of  m  values,  defined  in  (40), 
as  an  argument  for  the  elliptic  integrals  in  the  next  line.  The 
function  ZERO  replaces  each  number  whose  magnitude  is  less  than 
IE' 14  by  zero. 

[3]  K«-(AE1  MPOLY  M)-(©M)xBEl  MPOLY  M  O  E<-(AE2  MPOLY  M) 

-(@M)xBE2  MPOLY  M 

MPOLY  is  a  function  whose  left  argument  is  a  vector  of  polynomial 
coefficients.  It  evaluates  the  polynomial  for  each  element  in  M 
and  returns  the  results  in  the  same  shape  as  M.  AEl,  BEl,  AE2, 
and  BE2  are  the  coefficients  in  the  polynomial-and-log  approxima¬ 
tions  for  the  elliptic  integrals  given  in  [7  p591-592,  5  terms]. 
©M  gives  the  natural  log  of  the  values  in  M. 

[4]  G0[OX;SX]*-2xKvSQAB^AB=t<0.5 

O  G1  [0X;SX]<-  (2vB)  x  (Kx  (AB-B)  vSQAB)  -ExSQAB 

Line  [4]  expresses  and  assigns  the  blocks  of  element  values  to 
the  indexed  parts  of  GO  and  Gl.  The  first  statement  is  a  direct 
copy  of  (25).  The  second  statement  is  a  reworking  of  (28)  to 
minimize  the  number  of  variables  holding  intermediate  data. 

REALGREEN  called  its  geometry  values  by  indexing  the  previ¬ 
ously-constructed  vectors  for  p  and  z .  The  average  value  for 
each  G  function  over  a  singular  interval  requires  that  t'  values 
be  calculated  within  the  interval  for  the  integration  formula, 
and  the  G  functions  evaluated  at  these  t'  values.  Therefore  one 
can't  use  REALGREEN  and  some  of  the  same  statements  must  be 
rewritten.  This  is  done  in  the  function  IG. 

[0]  Z^IG  Y;X;K;E;A;B;AB;SQAB;XQ;M 

Z  and  Y  are  output  and  input  arguments  for  the  function,  each  a 
two-element  vector  as  can  be  seen  from  its  use  in  line  [6]  of 
SGENERATE.  The  local  variables  listed  in  the  header  mostly  have 
the  same  uses  as  they  do  in  REALGREEN. 
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[1]  X«-GIXxDTv2  O  XQ^XXX  O  B^2xY [ l]xY [ l]+XxY [2 ] 

O  AB<-XQ+2xB  O  A<-XQ+B 

GIX  is  a  vector  holding  the  tabulated  x,-  values  from  [7  p92  0, 
n=4].  DT^-2  is  the  half-interval  scale  factor  to  convert  the 

X,-  to  t'-t  values.  The  statement  for  B  is  a  copy  of  (37), 
letting  execution  order  take  care  of  the  parentheses.  AB  is 
(38).  A  is  (36).  These  are  all  four-element  vectors. 

[2]  M*-XQtAB  O  K<-(AE1  MPOLY  M)-(©M)xBE1  MPOLY  M  O  SQAB<-AB*0.5 
M  is  (40)  again,  and  K  and  SQAB  are  as  in  REALGREEN. 

[3]  E<-(AE2  MPOLY  M)-(©M)xBE2  MPOLY  M  O  Z*-2x+/GIWxK^SQABxLGX 

GIW  is  a  vector  holding  values  and  LGX  is  a  vector  holding 
log(Xj)  values.  '+/'  does  the  summation  so  now  Z  holds  the 
average  of  Gqe  over  a  particular  interval  for  which  it  is 
singular. 

[  4  ]  Z«-Z  ,  2X+/GIWX  (  (AXKvSQAB)  -SQABXE)  ^BXLGX 

Line  [4]  calculates  the  average  for  G^^  and  catenates  it 
with  Z  to  form  the  new  Z  which  the  function  passes  out  as  the 
explicit  result.  The  28  in  (33)  doesn't  appear  in  this  function 
because  (33)  is  the  integral,  not  the  average.  The  average  is 
the  integral  divided  by  25. 

The  complex  G  functions,  Gq^^  and  G^^  (23,24),  are  initial¬ 
ized,  collected  and  timed  in  a  function  WGENERATE  which  is  nearly 
identical  to  SGENERATE.  The  results  are  held  in  arrays  GOC  and 
GIC.  The  actual  integrations  are  done  by  the  function 
CMPLXGREEN. 

[0]  OX  CMPLXGREEN  SX;N;DA;ANG;R;CR;SR 

[1]  N<-15  O  DA*-0vN  O  ANG<-2oDAx‘0.5+iN 

O  R<-(RHT[OX]xRHT[OX])®  .+RHT[SX]xrhT[SX] 

The  integration  in  0  is  done  by  the  piecewise-constant  approxi¬ 
mation.  N  is  the  number  of  intervals  between  0  and  tt  ,  and  DA  is 
the  interval  length.  ANG  is  the  vector  of  cosine  values  for  the 
centers  of  these  intervals.  R,  in  this  line,  is  a  matrix  of  all 

the  possible  p^+p'^  values  specified  by  the  OX  and  SX  index 
vectors . 

[2]  R<-R+(ZT[0X]®  .-ZT[SX]  )  >^2 

<0  R<-(  (  (N,pR)  pR)-2xANG®  .xRHT[0X]®  .XRHT[SX])  *0.5 
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The  last  statement  in  line  [1]  and  these  two  statements  are  a 
three-step  construction  of  all  the  possible  values  of  the  dis¬ 
tance  function  required  by  the  index  vectors  OX  and  SX  and  the 

steps  in  0.  Coming  into  line  [2]  R  has  all  the  p^+p'^  values, 

and  the  first  statement  adds  all  the  (z-z')^  values.  Starting 
at  the  right  of  the  second  statement  in  line  [2],  the  ()*0.5  will 
take  the  square  root  of  all  the  elements  sent  to  the  outer 
parentheses.  Just  inside  the  right  parenthesis  pp'  is  done  as 
an  outer  product  to  form  a  matrix  the  same  shape  as  the  current 
value  of  R.  This  matrix  is  then  outer-product  multiplied  by  the 
COS0  values  in  ANG  to  form  a  3-axis  array.  The  first  axis 
corresponds  to  the  steps  in  <p,  the  other  two  to  the  index  values 
in  OX  and  SX.  This  array  is  next  multiplied  by  2  to  give  the 
2pp'cos0  part  of  the  distance  function.  Now  a  copy  of  the 
present  value  of  R  needs  to  be  generated  for  each  step  in  0  so 
the  subtraction  part  can  be  properly  done.  (N,pR)pR  makes  a 

3— axis  array  out  of  the  p^+p^^+(z— values  with  N  identical 
copies  along  the  first  axis. 

[3]  GOC[  ;OX;SX]«-DAx(+/CR<-('l+2oBTxR)^-R)  ,  [ 0 . 5 ] +/SR<- ( loBTxR)  tR 

[4]  ANG<-^(  ((ppR)  pANG)  O  GIC  [  :  OX;SX]  <-DAx  (+/ANGxCR)  ,  [  0 . 5  ] +/ANGxSR 

is  represented  in  rectangular  form  as  cosO^)  (2oBTxR)  and 
-sin()SR)  (-loBTxR)  .  The  imaginary  and  real  parts  are  separately 
formed  and  summed  down  the  first  axis  (the  0  steps)  and  lastly 
multiplied  by  DA  to  complete  the  integration.  The  first 
statement  in  line  [4]  reshapes  the  cos0  vector  into  a  3 -axis 
array  to  match  the  other  items  in  the  integrand  of  G^^^. 

The  reader  may  feel  that  these  are  compact  (a  good  thing) 
and  formidable  bits  of  code.  They  are  to  the  developer  too. 

They  were  developed  by  using  vectors  with  a  few  elements  each  to 
test  the  operations  and  see  that  they  worked  out  to  the  shapes 
needed.  If  an  operation  works  for  little  vectors,  it  will  work 
for  any. 

IVD  Filling  the  Impedance  Matrices 


Following  Simpson^s  path,  the  Z,-,^  matrix  was  split 
into  circuit-element  and  radiation  parts. 


2ik  ~ 


(47) 
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The  L  and  S  matrices  depend  strictly  on  geometry  with  no  frequen¬ 
cy  component,  so  that  once  they  are  computed  for  a  given  shape, 
they  need  not  be  recomputed.  It  actually  costs  more  execution 
time  to  compute  Z  this  way  most  of  the  time,  but  it  was  interest¬ 
ing  to  see  the  relative  effects  of  the  static  and  dynamic  G 
functions. 


L,k  = 


Sik  = 


f  i  Tk  ( Goecosecose '  +G,  E  s  inOs  inG ' )  dt '  dt 

Tj-A  t,j-A 


-J  J 


=  f  J 


r..6 


T,T|(  (GowCos0cos0'+Gi(jSin0sin0' ) 


(48) 

(49) 


dT;  Gq^  dT„ 
dt 


Idt'dt 


(50) 


Again,  two  kinds  of  functions  were  written  for  these  matrices, 
one  to  initialize  them,  set  up  loops,  and  time  the  calculations, 
and  one  to  do  the  actual  double  integration  given  the  i,k  values. 
[0]  REACT;I;K;Tl;NC 

This  function  sets  up  the  reactance  matrices,  L  and  S.  They  are 
called  XGL  and  XGC  respectively.  The  function  also  finds  the 
terminal  (circuit-element)  reactance  due  to  the  static  poten¬ 
tials. 


[1]  Tl^-^TS  O  1^1  O  XG«-XGL«"XGC^- (NP,NC*-1+2^NP)  PO 

[2]  IL:K^1 

These  lines  initialize  the  starting  time,  loop  counters,  and 
matrices.  NC  is  the  'N'  summation  limit  in  (46) ,  which  is  the 
number  of  equations  and  unknowns  without  taking  account  of  the 
geometric  symmetry. 

[3]  KL:XGL[I;K]^LIK  I,K  <>  XGC  [  I ;  K]  ♦‘SIK  I,K 
^  -»(NC^K^K+1)/KL  ^  -»(NP^I^I+1)/IL 

LIK  and  SIK  are  the  functions  that  do  the  integrations. 

[4]  XGL^(  (NP,NP)^XGL)+^0,  (NP,  1-NP) ^XGL 

[5]  XGC"-(  (NP,NP)^XGC)+<I^0,  (NP,  1-NP)  ^XGC 
O  XG^"XGL+XGC^BS 
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Here  is  where  the  symmetry  in  t  is  used  to  reduce  the  matrices  to 
NP  by  NP.  Suppose  NP=5.  Then  2  J9=Ji.  This 

implies  that  the  columns  to  the  right  of  5  can  be  reversed  in 
order  and  added  to  those  to  the  left  of  5.  Only  NP  rows  are 
calculated,  with  this  reduction  in  mind.  NP,1-NP  specifies  NP 
rows  and  NP-1  columns  from  the  right  side  of  the  matrix  for  the 
take  (T)  function.  The  0,  on  the  left  will  put  a  column  of 
zeros  on  that  side  of  the  NP  by  NP-1  present  result.  The  reverse 
flips  the  matrix  side-to-side  so  the  column  of  zeros  lines  up 
with  column  NP  and  column  NP+1  of  the  original  matrix  lines  up 
with  column  NP-1.  The  last  step  in  the  reduction  process  is 
addition  to  the  left  NP  columns  of  the  original  matrix.  XG  is 
the  total  reactance  matrix. 


[6]  JT< — (  (NP-1)  pO)  ,v60xBTxDTxDT)E1XG  O  XTERM«-vo2xJT [NP] 


JT  is  the  vector  of  current  density  coefficients,  J^|^.  The 

are  NP-1  zeros,  and  a  negative  imaginary  number,  -Iv  (eOySDt^)  . 
The  Dt  was  put  with  V,-  to  reduce  the  multiplies.  Since  V,-  is 
imaginary  and  XG  is  real,  JT  is  imaginary,  as  it  should  be. 
Vj,=l,  so  the  terminal  reactance  is 


^term 


27rp(-f) 


t„.T„p(^) 


NP  V  2 

A  /  \ 


'tNP 


1 

27rJt(jp 


(51) 


The  next  two  lines  of  the  function  print  the  reactance,  the  run 
time,  and  sound  the  scale. 

The  integration  functions  have  the  same  format,  so  only  the 
more  interesting  LIK  is  described. 

[0]  Z*-LIK  IK;L;M 

[1]  L*-P+2^IK[1]  O  M«-P+IK[2] 

[2]  Z«-(  (TRT^SINTCL]  )  +  .><Gl[L;M]  +  .^SINT[M]^TRT)  + 

(TrtXcoST [ L] ) + . ^GO [ L;M] + . ^COST [M] ^TRT 

P  is  a  global  vector  whose  elements  are  '1  0  1  2.  This  makes  L 
and  M  the  set  of  t  and  t'  indices  needed  for  the  four-point  inte¬ 
gration,  as  shown  in  Figure  4.  The  appropriate  portions  of  the 
GO  and  G1  matrices  are  called  as  4x4  matrices.  Each  one  is  then 
the  center  matrix  in  a  vector-matrix-vector  multiply  (inner 
product)  to  do  the  double  sums.  That's  it. 
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Another  pair  of  functions,  IMP  and  ZIK,  organize  and  compute 
the  elements  for  The  computation  principle  for  the  integration 

is  the  same  as  in  LIK,  but  the  complex  numbers  require  a  little 
more  work  in  handling  the  shapes  of  the  vectors  and  arrays.  The 
reactance  and  impedance  matrices  are  combined  and  the  net 
terminal  impedance  computed  in  the  function  TERM. 

IVE  Current  and  Dissipation 

Once  the  current  density  coefficients,  Jj,^,  are  known, 
the  current  at  each  sample  point  can  be  found  by  adding  the 
appropriately-weighted  triangle  functions.  The  values  of 
are  held  in  the  variable  JTC,  but  only  for  half  the 
dipole.  To  provide  plots  for  the  exercise  of  intuition,  the 
current  was  extended  for  the  whole  dipole.  All  of  the  A  seg¬ 
ments  have  two  triangles  on  them  except  the  starting  one  at  t=0 
and  the  ending  one  at  t=t^.  The  total  current  is  formed  by 
generating  a  complex  vector  of  coefficient  values  padded  with  0 
at  the  right  end  (t^)  to  multiply  by  the  first  half  of  the  tri¬ 
angle  function,  and  another  coefficient  vector  padded  with  0  at 
the  left  end  {t=0)  to  multiply  by  the  second  half  of  the  triangle 
function,  and  then  add  the  two. 

ICX*-(JTC,(0  0  ‘1  iJTC),0  0)®.xTRT[l  2]  (52) 

Since  JTC  is  a  complex  vector  (2-row  matrix)  the  outer  product 
produces  a  two-plate,  or  3-axis,  result.  The  first  column  of  the 
first  plate  is  0.25  times  the  real  parts,  the  second  column  of 
the  first  plate  is  0.75  times  the  real  parts,  and  the  second 
plate  is  the  corresponding  imaginary  parts  of  the  current  coeffi¬ 
cient  times  0.25  and  0.75. 

RC*-(,ICX[1;;]) ,  [0.5]  ,ICX[2; ;  ]  (53) 

After  extending  JTC,  there  are  2NP  by  2  elements  in  each  plate  of 
ICX.  Raveling  each  plate  takes  the  plate  row-by-row  and  makes 
one  4NP-element  row,  every  other  element  is  0.25  times  a  coeffi¬ 
cient  and  its  neighbor  to  the  right  is  0.75  times  that  same 
coefficient.  Thus,  each  pair  in  a  row  is  the  left  half  of  a 
triangle  function  times  a  current  coefficient.  The  next  state¬ 
ment  reuses  ICX  to  multiply  the  current  coefficients  by  the  right 
half  of  the  triangle  function,  with  an  extra  zero  put  in  at  the 
beginning.  In  (55)  ICX  is  reformatted  as  above  and  the  two  sets 
are  added. 

ICX«-(0  0  ,JTC,0  0  ’liJTC)  o  .xTRT[3  4]  (54) 

RC^RC+(,ICX[1; ;]) , [0.5] ,ICX[2; ;]  (55) 
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RC  holds  the  surface-current  density  times  the  radius,  because 
the  complete  basis  function  hasn't  been  used.  To  get  the  total 
current  through  a  cross-section,  just  multiply  RC  by  Iw . 

IT<-02XRC  (56) 


One  of  the  banes  of  electrically-small  antennas  is  that  loss 
resistance  can  be  higher  than  radiation  resistance.  One  of  the 
benefits  expected  from  a  thick  dipole  is  a  lower  loss  resistance 
due  to  spreading  the  current  out  over  a  large  surface.  Assuming 
that  which  is  desired,  the  dissipation  should  have  little  effect 
on  the  current  distribution,  so  one  can  use  the  current  from  the 
lossless  equation  to  calculate  the  dissipation.  The  total 
dissipated  power  is 


P  =  27rJ^'"RsJtJtPdt 
with  Rg  = 


^  _  te60F 

“  \  a 


2a 


(57) 

(58) 


At  this  point  in  the  programming,  RC  is  J^P,  so  magnitude  squared 
of  RC  divided  by  RHT  (P)  is  used  instead  of 


p«-02XDTX  (  (oeOXBT^SIGMA)  *0.5)x+/  (+/RCXRC)  ^RHT  (59) 


V  Verification 

One  cannot  build  a  system  of  APL  functions,  any  more  than 
one  can  do  a  long  mathematical  development,  test  it  only  at  the 
end,  and  hope  to  live  without  ulcers.  Verification  begins  by 
testing  the  operation  of  each  function  on  simple  cases  that  can 
be  verified  by  inspection,  by  hand  calculation,  or  by  looking  up 
results  in  tables.  Also,  limiting  values  for  large  or  small 
arguments  are  sometimes  useful  checks.  Both  the  Green's  function 
matrices  and  the  impedance  matrices  (before  folding)  must  be 
symmetrical.  In  earlier  stages  of  this  project,  more  of  the  work 
was  done  by  completely  separate  functions,  and  there  were  more 
loops  in  the  matrix-generating  lines.  This  allowed  easier 
testing  of  things  like  the  elliptic  integrals  and  their  argument. 
Once  some  confidence  was  developed  in  results  from  these  simpler 
functions,  these  results  were  kept  for  comparison,  in  separate 
workspaces,  against  results  from  more  sophisticated  versions. 

The  final  tests  are  overall  performance  tests.  Since  we  are 
approximating  a  continuous  current  distribution  by  a  finite  set 
of  sample  values,  and  using  approximate  integrations  everywhere, 
we  cannot  expect  perfect  accuracy.  We  have  a  right  to  expect 
that  the  terminal  impedance  values  will  behave  smoothly  as  a 
function  of  the  number  of  sample  points,  and  converge  to  some 
result.  The  slow  convergence  of  the  impedance  values  was,  at 
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first,  daunting.  Reading  about  the  similar  experience  of  others 
[14,  15]  gave  a  little  morale  boost.  Although  guite  a  number  of 
test  shapes  were  tried,  with  varying  success,  two  are  offered 
because  they  are  shape  extremes  and  are  represented  in  the 
literature.  These  are  the  very  thin  wire  dipole  and  the  sphere 
dipole,  both  operated  at  half  wavelength. 

The  classical  impedance  for  a  thin  halfwave  dipole  is 
73+j42.5  n.  To  describe  this  shape  to  the  functions,  set 
RHD=1E'8  0,  ZD=0.9999  1.  There  was  no  difference  in  the  results 
for  the  radius  two  orders  of  magnitude  larger.  The  triangle 
basis  function  automatically  sets  the  current  at  the  dipole  end 
to  zero,  and  the  difference  between  the  two  ZD  values  is  much 
less  than  any  integration  segment  used,  so  the  result  is  a 
numerical  model  for  a  very  thin  straight  wire.  The  ZD  values 
couldn't  be  identical  because  a  divide-by-zero  error  would  occur. 
Because  of  experience  in  control  systems  and  broadband  rf  elec¬ 
tronics,  convergence  was  originally  tested  by  doubling  NP  from 
one  trial  to  the  next,  later  switching  to  a  2-5-10  sequence.  The 
largest  number  which  could  be  run  without  overfilling  the  work¬ 
space  was  NP=150.  This  corresponds  to  600  sample  points. 

Finally,  NP=5  10  20  50  100  was  settled  on  as  the  test  values.  The 
admittance  values  were  then  fitted  by  polynomials  in  1/NP.  This 
course  was  taken  because  a  polynomial  in  1/NP  goes  to  a  constant 
value  as  NP-»<».  Thus,  if  the  solution  values  do  converge,  the 
polynomial-fit  constants  might  be  the  right  values.  Table  1 
shows  the  data  and  results  for  the  thin  dipole. 

Table  1.  Admittances  for  the  Thin  Halfwave  Dipole,  mS 

NP  5  10  20  50  100 

G  11.05  10.13  9.823  9.709  9.7 

B  -4.548  -5.071  -5.278  -5.42  -5.486 

G  =  9.713-2. 501NP'''+133NP'^-889.2NP'^+2270NP’^ 

B  =  -5.568+9.077NP’''-99.66NP'^+772NP'^-1866NP'^ 

Z„  =  77.49  +  j44.42  D 

The  spherical  dipole  profile  has  to  be  approximated  by  a 
sequence  of  straight-line  segments.  A  five-segment  model  with 
equal  subtended  angles  was  chosen  for  the  half-sphere  so  that  the 
longest  A  value  would  fit  between  the  corners.  The  generating 
statements  are  ZD<-0,  looO .  Ixi5  and  RHD<-l,2ooO.  Ixi5.  The  results 
for  the  same  trial  values  and  curve  fitting  are  given  in  Table  2. 


00 

9.713 

-5.568 
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Table  2.  Admittances  for  the  Halfwave  Sphere  Dipole^  mS 


NP 

5 

10 

20 

50 

100 

00 

G 

17.01 

16.71 

16.65 

16.64 

16.64 

16.64 

B 

58.68 

71.41 

83.78 

99.96 

112.2 

130.1 

G  =  16.64-0.1826NP'’+5.024NP*^+47.95NP’^-112.1NP‘^ 


B  =  130.1--2119NP'''+3.592E4NP’^-2.764E5NP'^+7.043E5NP‘^ 

Z«  =  0.9673  -  j7.563  n 

The  polynomial  coefficients  give  an  indication  of  the  rate 
of  convergence,  if  it  exists.  The  sphere  represents  two  ex¬ 
tremes,  a  case  where  convergence  has  already  occurred,  and  a  case 
for  which  it  is  doubtful.  NP  has  to  go  over  2000  to  drive  the 
variable  terms  below  1  mS  for  B.  For  the  thin  dipole,  NP=200 
drives  the  variable  part  of  B  below  1  %  of  the  constant,  and 
NP=25  will  do  the  same  for  G. 


Figure  6.  Cross-sectional  current  along  the  half-wave  thin 
dipole. 
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Figure  7.  Cross-sectional  current  along  the  half-wave  sphere. 

Forging  ahead,  the  current  distributions  were  plotted  for 
spheroids  of  various  radii  froin  1  m  on  down.  Figures  6  and  7 
show  the  two  extremes.  The  real  parts  show  a  sinusoidal  shape  as 
might  be  expected.  The  imaginary  part  of  the  thin  dipole  current 
is  negative,  necessary  for  the  inductive  terminal  impedance,  and 
looks  sinusoidal  except  for  a  tiny  upward  cusp  at  the  feedpoint. 
As  the  dipole  is  thickened,  this  cusp  grows,  the  imaginary  part 
of  the  current  changes  sign  along  the  z  axis,  and  finally  becomes 
entirely  positive.  The  full  sphere  shows  a  large  leading  current 
near  the  feed  gap.  Since  this  current  is  a  growing  function  of 
the  antenna  radius,  it  is  presumed  that  it  is  bounded  as  it 
looks,  and  the  admittance  values  are  reasonable,  if  not  highly 
accurate. 

According  to  Ramo  and  Whinnery,  Stratton  and  Chu[ll]  used 
spherical-harmonic  modes  to  match  the  boundary  conditions  for 
spheroidal  antennas.  Their  results  for  the  full  sphere  give  a 
converging  series  for  the  conductance,  but  not  for  the  suscep- 
tance.  Figure  12.23d,  p544,  from  [12]  shows  curves  for  these 
quantities  as  functions  of  ySa,  'a'  being  the  sphere  radius.  The 
susceptance  value  for  )3a=7r/2  from  this  figure  is  about  27  mS, 
which  was  expected  to  be  low,  because  they  stopped  arbitrarily  at 
19  terms.  The  conductance  value  is  about  17  mS,  which  matches 
Table  2  above.  A  frequency  response  for  )8a  between  0  and  3  was 
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VI .  Impedance  and  Shape 

For  many  years  it  seemed,  from  [12],  that  the  sphere  was  a 
potential  broadband  antenna,  but  the  recent  reexamination  of  this 
material  and  the  calculations  described  in  section  V  show  that 
the  large  gap  perimeter  adds  too  much  capacitance  to  the  antenna. 
At  small  a/\  the  wire  dipole  has  too  much  series  body  capaci¬ 
tance  and  the  sphere  has  too  much  feed  capacitance.  A  shape  with 
more  body  than  the  wire  and  less  gap  than  the  sphere  ought  to  do 
better  than  both.  From  Hansen's  paper[3]  quoting  Chu[13],  the 
minimum  Q  for  an  antenna  with  /3a<l  is 


=  l+3i6^a^ 


(60) 


At  3  MHz,  K=100  m,  and  the  antenna  size  limit  is  a=l  m.  This 
means  that  Q^,.^=4063.  For  a  simple  tuned  circuit  this  means  a 
bandwidth  of  only  about  750  Hz,  a  little  tight  even  for  SSB. 

This  may  seem  a  high  value  of  Q  but  no  reported  wire  antennas 
come  anywhere  near  this  low.  Starting  from  the  sphere,  the  gap 
was  opened  up  by  using  a  cone  from  a  small  gap  as  a  transition 
section  to  a  barrel  and  a  capping  cone,  without  good  result. 

Then  a  biconical  half-dipole  was  explored  and  showed  that  the  Q 
improved  as  the  cap  cone  was  flattened.  This  led  to  a  single 
cone  flaring  out  from  the  feed  to  a  plate  cap.  Is  the  straight- 
sided  cone  the  best  that  can  be  done?  Both  convex  and  concave 
profiles  were  tried.  Concave  was  better.  Finally,  the  tube  and 
cone  was  found  to  be  the  best  of  this  shape  type.  The  profile 
for  the  best  tube-and-cone  half-dipole  is  a  tube  of  10  cm  radius 
from  the  feed  gap  to  0.75  m,  then  a  cone  flaring  out  to  1  m 
radius  at  the  end,  with  a  plate  to  close  the  figure.  Table  3 
gives  descriptions  and  results  for  some  of  the  trial  shapes. 


A  basic  conclusion  one  can  draw  from  this  work  is  that  both 
gap  and  body  capacitances  need  to  be  kept  down.  Even  with  the 
gap  perimeter  reduced  going  from  the  sphere  to  the  barrel, 
there's  a  great  deal  of  charging  current.  So  why  isn't  a  tube 
out  to  a  plate  the  best  shape?  Possibly  a  little  circuit  tuning 
is  going  on.  A  frequency-response  run  for  the  T-C  (tube  and  cone 
antenna),  shown  in  Table  4,  turned  up  a  series  resonance  near 
17.49  MHz. 


Figures  9-12  show  the  current  distributions  along  the  T-C 
dipole.  The  real  part  is  nearly  flat  along  the  tube  below,  at, 
and  above  resonance,  while  the  imaginary  part  is  closer  to  flat 
than  anything  else  below  and  above  the  resonance.  This  is  what 
one  would  expect  from  an  end-loaded  small  dipole  and  supports  the 
idea  that  the  resonance  is  a  circuit  phenomenon.  The  sharp  drop 
in  the  z  components  near  the  dipole  ends,  shown  in  Figure  10,  are 
due  to  the  abrupt  change  in  the  current's  direction  going  from 
the  tube  to  the  wide-angle  cone.  The  current  shifts  from  axial 
to  almost  radially  directed. 
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Table  3 .  Circuit  Properties  for  Dipole  Shapes 


f=3  MHz,  N=100  m,  dipole  length=2  m,  NP=20.  and  are 
the  half-dipole  profile  description  in  meters.  Rioss 
Z,-^  are  series  equivalent  circuit  parameters.  Q=|X|/R. 

Shape  Name  pyz^  Z,-„  n  Q  Ri^ss  mfl  Cu 


Thin  Spheroid 
2  mm  feed  dia. 

Full  Sphere 

Barrel  0.01  1  1  0 

0  0.2  0.9  1 

Double  0.01  1  0 

Cone  0  0.5  1 

Cone  and  0.01  1  0 

Cap  0  0.99  1 

"  0.01  0.5  1  0 

0  0.5  0.999  1 

Concave  0.01  0.5  1  0 

Cone  0  0.8  0.999  1 

Tube  and  0.01  0.01  1  0 

Cone  0  0.75  0.999  1 

(T-C) 


0.067-jll,290 

168k 

44.8 

0.0142-j267 

19k 

0.0128 

0.019-j233 

12.1k 

0.435 

0.043-j504 

11.7k 

0.513 

0. 142-j673 

4765 

0.672 

0. 145-j679 

4685 

0.677 

0.188-j807 

4293 

0.895 

0.269-j925 

3432 

10.4 

Table  4 .  Frequency  Response  of  the  T-C  Dipole 


f,  MHz 

3 

6 

10 

Zin  n 

0.27-j925 

1.09-j422 

3.1-jl95 

f,  MHz 

17.49 

20 

30 

a 

c 

10.28+j0.073 

14+j47 

39+j215 
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run,  and  Figure  8  shows  the  plot  of  conductance  for  /fe  between  0 
and  2.5.  The  parameters  were  NP=20  and  the  half-sphere  profile 
was  divided  into  20  equal  segments.  The  curve  has  the  same  shape 
as  that  from  [12].  That's  nice,  but  why  aren't  all  the  conduc¬ 
tance  data  plotted?  Because  there's  a  big  wrinkle  around  /fe=2.8 
which  would  have  obscured  the  shape  of  the  lower-frequency  curve. 
The  sphere  has  a  series  resonance  at  about  ^=2.752,  K=2.283, 
and  a  parallel  resonance  at  about  )9a=2.882,  K=2.18  m.  Tinkering 
with  NP  and  segment  number  did  not  make  the  resonances  go  away  so 
they  appear  to  be  real.  Fishing  for  the  series  resonance  K  was 
interesting  because  the  susceptance  behaved  reasonably,  but  the 
conductance  did  strange  things  in  magnitude  and  changed  signs 
close  to  the  resonance.  Since  the  conductance  should  get  large 
and  the  susceptance  small  at  series  resonance,  perhaps  the 
impedance  matrix  was  becoming  difficult  to  invert.  The  parallel 
resonance  is  smoother. 


Figure  8.  Conductance  for  a  40-segment  sphere  dipole. 

Since  agreement  with  two  classical  cases  of  very  different 
shape  is  fair,  it  seems  reasonable  to  gamble  that  the  functions 
and  numerical  procedures  are  generally  working  and  reliable. 

Ramo  and  Whinnery[12]  said  the  series  for  the  susceptance  doesn't 
converge  for  the  infinitesimal  gap,  so  the  mode  analysis  has  the 
basic  defect  that  it  can't  account  for  the  circuit-element  behav¬ 
ior  of  the  sphere.  This  may  be  why  the  mode  analysis  doesn't 
show  up  the  resonances. 
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Figure  11.  Current  on  the  T-C  dipole  at  17.49  MHz.  Real  part  in 
A,  imaginary  part  in  dA. 


Figure  12.  Current  on  the  T-C  dipole  at  30  MHz. 
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VII.  Development  Time  vs.  Execution  Time 


This  project  was  spread  over  a  year  and  a  half  beginning  in 
January,  1989.  Too  much  time  was  spent  in  the  first  quarter  of 
1990  in  sorting  out  hardware  and  software  problems  in  a  move  to  a 
32-bit  pc  and  a  32-bit  APL  interpreter.  Writing  this  paper  has 
taken  three  7 -day  weeks.  By  comparison,  the  time  spent  in  the 
code  development  seems  only  a  brief  pleasant  memory.  The  reader 
has  seen  almost  all  the  code  written  for  this  project  in  this 
paper.  Think  about  that.  Realize  how  closely  the  code  mirrors 
the  mathematics  and  how  little  overhead  there  is.  In  1989,  an  8 
MHz  286/287  system  with  2.7  Mb  of  RAM  was  being  used  for  this 
project.  A  few  times  cases  were  run  that  took  about  5  hours. 

With  the  move  to  the  20  MHz,  8  Mb,  32-bit  system,  the  longest 
case  took  about  an  hour.  The  100-point  admittance  calculation 
for  Figure  8  took  about  2  hours,  during  which  markup  on  a  stu¬ 
dent's  thesis  draft  was  started.  Now  this  project,  which  is  part 
of  a  larger  one,  is  over.  This  code  may  never  run  again.  From 
this  perspective,  even  if  the  execution  time  of  a  compiled 
program  was  zero,  it  would  not  compensate  for  the  excess  time 
spent  in  coding  and  debugging  the  constituent  functions  in  a 
scalar  programming  language.  Human  time  is  too  precious. 

Table  5  shows  the  execution  times  for  various  parts  of  a 
case  computation.  The  blank  entries  are  for  cases  that  won't  fit 
in  8  Mb.  While  these  times  are  for  a  particular  case,  a  simple 
shape  and  long  wavelength,  there  is  considerable  variation  from 
case  to  case.  Shorter  wavelengths  seem  to  take  longer,  as  much 
as  10  %  on  the  dynamic  calculations.  This  effect  has  not  been 
investigated.  The  array  sizes  don't  depend  on  wavelength,  nor 
are  there  any  loops  which  test  on  relative  size  for  termination. 
The  maximum  problem  that  will  fit  with  row  looping  is  NP=150  (600 
samples) ,  and  looping  on  ten-element  bites  out  of  a  row  allows 
the  problem  to  run  on  a  640  kb  machine  with  the  16-bit  interpret¬ 
er.  The  times  for  the  G  functions  with  ten-element-bite  loops 
are  double  those  for  no  loops  on  the  same  machine.  One  can  see 
that  the  times  generally  follow  an  NP^  law.  The  difference 
between  total  case  time  and  time  for  a  new  frequency  is  a  little 
deceiving,  if  one  is  considering  whether  to  fill  a  separate 
static  matrix  or  not.  If  the  static  and  dynamic  G  functions  are 
added  and  then  used  to  fill  only  one  impedance  matrix,  the  time 
to  fill  that  matrix  would  be  only  slightly  longer  than  the  times 
given  for  the  dynamic  matrix  fill.  In  such  a  scheme,  the  time 
for  a  new  frequency  point  in  a  response  run  would  be  the  same, 
and  the  saving  per  frequency  would  be  the  time  to  compute  the 
static  G  functions.  The  circuit-element  reactance  was  calculated 
separately,  for  cases  with  new  geometry  or  number  of  samples,  for 
interest.  The  utility  of  this  data  is  unclear,  so  it  wasn't 
included  in  this  paper. 
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Table  5.  Execution  Times 

All  times  are  in  seconds.  Cases  were  run  on  an  ATronics  386B/20 
(20  MHz  clock)  w.  64k  cache  card,  8  Mb  RAM,  Intel  '387.  The 

release  2 . 


Sample  point,  4xNP. 

Static  Gq  and  without  loops. 
Static  Go  and  G^  with  row  loop. 
Static  matrix  fill  and  invert. 
Dynamic  Gq  and  G^  without  loops. 

Dynamic  Gq  and  G^  with  row  loop. 
Dynamic  matrix  fill  and  invert. 
Total,  without  G  loops. 

Total,  with  G  loops. 

Time  for  a  new  frequency. 


20 

40 

80 

200 

400 

1.16 

3.4 

11.3 

63.5 

1.49 

4.12 

13.5 

66.8 

249 

1 

4 

16.2 

103 

426 

2.9 

11.4 

45 

3.1 

12 

47.4 

284 

1140 

1.4 

5.1 

22.2 

157 

789 

6.5 

24 

95 

7 

25 

100 

511 

2704 

in 

• 

17 

70 

441 

1929 
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Appendix;  APL  Syntax,  Symbols  and  Functions 

This  appendix  will  present  only  as  much  APL  as  is  used  in 
the  paper.  APL  is  fundamentally  an  interactive  array-processing 
language.  To  be  interactive,  it  must  be  an  interpreter.  It  is 
economical  of  time  in  development  and  execution  because  it  uses 
an  extended  symbol  set  from  which  over  80  functions  and  operators 
can  be  called  by  one  or  two  keystrokes  each.  There  are  no  data 
type  statements,  a  variable  is  given  its  type  by  the  way  the  user 
assigns  data  to  it,  and  the  interpreter  keeps  track  of  the  type, 
not  the  user.  Arithmetic  is  double-precision,  or  better.  A 
variable  may  be  a  scalar,  a  vector,  or  an  array  of  any  number  of 
axes  (dimensions) .  Execution  of  a  statement  proceeds  from  right 
to  left,  except  when  interrupted  by  parentheses.  AxB+C  is  not 
equal  to  C+AxB.  Thus,  a  typical  statement  execution  begins 
with  the  interpreter  reading  data  (explicitly  or  in  a  variable 
name)  at  the  right  end  and  finishing  either  with  a  screen  display 
or  an  assignment  to  a  variable  at  the  left  end.  Customarily, 
only  uppercase  letters  are  used  in  names  for  variables  and 
user-defined  functions,  but  the  interpreter  is  case-sensitive  so 
that  lowercase  letters  can  also  be  used.  The  O  is  a  statement 
separator,  so  that  more  than  one  statement  can  be  placed  on  the 
same  line.  The  intepreter  passes  from  one  statement  to  the  next 
in  the  usual  left-to-right  manner,  as  if  the  succesive  statements 
did  have  line  numbers.  Multiple-statement  lines  is  a  feature  of 
STSC's  APL^PLUS  systems,  but  it  is  not  standard.  In  the 
following  listing,  related  symbols  are  grouped  together  to  save 
space.  S  is  a  scalar,  N  is  an  integer,  V  is  a  vector,  M  is  a 
matrix,  f  and  g  represent  any  built-in  (primitive)  function. 
Primitive  scalar  functions  operate  on  a  variable  in  an 
element-by-element  manner,  fM  does  f  to  each  element  in  M  and  the 
result  has  the  same  shape  as  M.  There  are  also  array  functions 
which  rearrange  arrays  or  extract  blocks  of  data  from  them. 
Operators  modify  the  way  functions  are  applied  to  data. 
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APL  Symbols 


Meanings 


+  -  X 


Basic  arithmetic  functions. 


*X  X*Y  ©X 
X*-Y 
V 
iN 


LABEL: 


e\  X\  ln(X). 

X  is  assigned  the  shape  and  values  of  Y. 

Open  or  close  (toggle)  function-definition  mode. 

A  vector  of  integers  from  1  to  N. 

Go  to.  -4(test)/line_label.  If  'test'  is  true 
control  passes  to  the  statement  line  beginning 
with  'line_label' . 

A  line  label  is  identified  by  the  colon. 


oX 

NoX 

V*-pX  A<-VpX 


|X  XX 

VEIM  |±1M 
,X  VI, V2 
VI, [0.5]V2 
ia  0  0 

VTM  ViM 

f/M  f/M 

Mlf .gM2 
XI® . fX2 


TT  times  X. 

N=1  is  sinX,  =2  is  cosX,  etc. 

Shape  is  the  number  of  elements  along  each  axis 
of  an  array.  pX  gives  the  shape  of  X,  a  vector. 
With  a  vector  left  argument  the  funtion  p  re¬ 
shapes  the  data  in  X,  taking  it  row  by  row,  to 
the  shape  specified  by  the  vector. 

Magnitude  of,  sign  of,  each  element  in  X. 

M'V,  M’''. 

Make  X  a  vector,  make  VI  and  V2  a  vector. 

Place  VI  over  V2  in  a  two-row  matrix. 

Matrix  rearrange.  Transpose,  reverse  left  to 
right,  reverse  top  to  bottom. 

Take  and  Drop.  V  specifies  how  many  rows  and 
columns  of  M  are  affected. 

Reduction.  Apply  f  between  the  elements  of  each 
row  (column)  of  M. 

Inner  product.  f/MlgiaM2. 

Outer  product,  f  is  applied  between  each  ele¬ 
ment  of  XI  and  every  element  of  X2.  The  shape 
equals  the  joining  of  the  shapes  of  XI  and  X2. 


131 


^0  f  037GI1.COS 

"[11  Mautz ,  J.  R.  and  R.  F.  Harrington,  "Radiation  and  Scattering 
from  Bodies  of  Revolution",  Applied  Scientific  Research  vol. 
20,  June  1969. 

[2]  Simpson,  T.  L,  J.  C.  Logan,  and  J.  W.  Rockway,  "Equivalent 
Circuits  for  Electrically  Small  Antennas  Using 
LS-Decomposition  with  the  Method  of  Moments", 

IEEE  Trans.  Ant,  and  Prop,  vol.  37,  no.  12,  December  1989. 

[3]  Hansen,  R.  C.,  "Fundamental  Limitations  in  Antennas", 

Proc.  IEEE  vol.  69,  no.  2,  February  1981. 

[4]  Miron,  D.  B.  "Survey  of  Recent  Results  for  the  Singularity 
Problem  in  the  Electromagnetic  Integral  Equation", 

Proc.  North  Dakota  Acad.  Sci.,  vol. 38,  p30,  April  1984. 

[5]  Mahadevan,  K.  and  H.  A.  Auda,  "Electromagnetic  Field  of  a 
Rectangular  Patch  of  Uniform  and  Linear  Distributions  of 
Current",  IEEE  Trans.  Ant,  and  Prop,  vol. 37,  no.  12, 

December  1989 

[6]  Gradshteyn,  I.  S.  and  I.  M.  Ryzhik, 

Table  of  Integrals.  Sums,  and  Products, 

Academic  Press,  New  York,  1965,  pl53-154. 

[7]  Abramowitz,  M.  and  I.  A.  Stegun, 

Handbook  of  Mathematical  Functions, 

Dover,  New  York,  1972(7),  p591-592. 

[8]  Pierce,  B.  0.  and  R.  M.  Foster,  A  Short  Table  of 
Integrals .  4th  ed. ,  p72,  Ginn,  New  York,  1956. 

[9]  Neff,  H.  P.  jr.,  Basic  Electromagnetic  Fields. 

2nd  ed..  Harper  and  Row,  1987,  pl54-159. 

[10]  Miron,  D.  B. ,  "The  Conducting  Cylinder  Problem",  an 
unpublished  handout,  available  from  the  author. 

[11]  Stratton,  J.  A.  and  L.  J.  Chu,  J.  AppI.  Physics,  no.  12, 
pp230-248,  March  1941. 

[12]  Ramo,  S.  and  J.  R.  Whinnery, 

Fields  and  Waves  in  Modern  Radio,  2nd  ed. , 

Wiley,  New  York,  1960. 

[13]  Chu,  L.  J.  "Physical  Limitations  of  Omnidirectional 
Antennas",  J .  AppI .  Phvs .  no.  19,  ppll63-1175,  Dec.  1948. 

[14]  Sarkar,  T.  K. ,  A.  R.  Djordjevic,  and  E.  Arvas,  "On  the 
Choice  of  Expansion  and  Weighting  Function  in  the  Numerical 
Solution  of  Operator  Equations",  IEEE  Trans.  Ant,  and  Prop. 
vol.  AP-33,  no.  9,  p988,  September  1985. 

[15]  Petersen,  A.  F.,  "Difficulties  Encountered  When  Attempting 
to  Validate  Thin-Wire  Formulations  for  Linear  Dipole 
Antennas",  Applied  Computational  Electromagnetics  Society 
Journal .  Special  Issue  on  Code  Validation,  p41,  1989. 


132 


ANALYSIS  OF  THREE  DIMENSIONAL 
DIELECTRIC  LOADED  CAVITIES 
WITH  EDGE  ELEMENTS 


L.  Pichon  A.  Razek 

Laboratoire  de  Genie  Electrique  de  Paris 
Ecole  Super ieure  d'Electricite 
U.R.A.  D0127  CNRS 
Universites  Paris  6  et  Paris  11 
91  192  Gif  sur  Yvette  Cedex 
France 


ABSTRACT 

In  this  paper  we  show  that  edge  elements  (a  class  of 
mixed  finite  elements)  provide  an  efficient  numerical 
approach  in  the  determination  of  resonant  modes  in  three 
dimensional  high  frequency  cavities.  These  finite  elements 
avoid  "spurious  modes”,  the  non-physical  numerical  fields 
obtained  from  the  solution  of  eigenvalue  problems. 

Here,  empty  cavities  as  well  as  dielectric  loaded 
cavities  are  analyzed;  no  "spurious  mode”  was  observed. 
Moreover,  comparisons  with  analytical  results  and 
previously  published  ones  show  the  great  accuracy  of  the 
numerical  technique. 

INTRODUCTION 

Electromagnetic  resonance  is  important  in  the  design 
of  particle  accelerators,  microwave  ovens  and  resonant 
cavities.  For  such  analysis,  numerical  techniques  including 
the  finite  element  method  have  been  developed. 

The  well  known  finite  element  method  seems  very 
attractive  since  for  several  years  it  has  been  found  to  be 
an  efficient  tool  in  low  frequency  electromagnetic  field 
computations.  In  high  frequency  applications,  finite 
elements  were  used  for  cavity  resonances  analysis  [l]-[5]; 
resonant  modes  and  resonant  frequencies  are  obtained  as 
solutions  from  an  eigenvalue  problem. 

The  main  serious  drawback  in  these  studies  is  that  the 
computed  solutions  are  plagued  by  non-physical  (  or 
"spurious”)  solutions:  solutions  which  do  not  satisfy  the 
divergence  free  condition  implied  by  the  Maxwell's 
equations.  Many  attempts  were  performed  to  circumvent  these 
unwanted  numerical  fields  (enforcing  the  divergence 
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condition  with  a  global  penalty  term  [1],[5],  reducing  the 
number  of  unknowns  by  taking  locally  into  account  the 
divergence  condition  [3]  ,  using  divergence  free  trial 

functions  ...  ) .  It  has  been  observed  that  discretized 
fields  with  continuous  tangential  components  suppress  the 
"spurious  modes"  (the  problem  was  solved  for  a  scalar 
function  in  two  dimensions  [5]  and  for  the  magnetic  vector 
potential  [4]  in  three  dimensions) ;  nevertheless  no  precise 
argument  was  put  forward  to  explain  the  importance  of  this 
kind  of  approximation. 

However  the  choice  of  finite  elements  for 
electromagnetic  field  computations  is  essential:  A. 
Bossavit  showed  that  "edge  elements"  are  well  adapted  for 
the  representation  of  vector  fields  since  they  allow  their 
possible  discontinuities  [6],  These  "edge  elements"  are 
also  a  class  of  the  mixed  finite  element  proposed  by  J.C. 
Nedelec  [7].  Such  elements  were  successfully  used  in  eddy 
currents  problems  [8] -[10]  and  are  well  adapted  for  the 
approximation  of  scattering  and  resonance  problems  [11] 
[12]  ;  the  reason  for  which  they  would  not  generate 

"spurious  modes"  is  explained  in  [12]. 

We  have  developed  and  applied  such  a  numerical 
approach  for  empty  and  dielectric  loaded  cavities.  In  this 
paper,  we  present  first  the  variational  formulation  of  the 
Maxwell's  equations  in  terms  of  electric  field  and  the 
reason  of  the  occurence  of  "spurious  modes".  Then  we  detail 
the  numerical  discretization  and  explain  the  interest  of 
"edge  elements".  Finally  we  present  the  analysis  of  three 
dimensional  cavities  . 

VARIATIONAL  FORMULATION 

We  deal  with  the  Maxwell's  time-harmonic  equations  in 
a  bounded  region  surrounded  by  a  perfect  conductor  and 
containing  lossless  materials  : 

rot  e  =  -  iti)  p.Qh  (1) 

rot  h=ia)€g€^e  (2) 

where  e  and  h  are  the  complex  electric  and  magnetic  fields, 
€q  and  are  respectively  the  permittivity  and 

permeability  of  vacuum,  is  the  relative  permittivity 

and  w  is  the  angular  frequency. 

The  conditions  on  the  boundary  r  of  fJ  are  those  of  a 
perfect  conductor: 

n-e  =  0  (3)  n.h  =  0  (4) 
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where  n  is  the  outward  normal  vector. 

Substituting  (1)  into  (2),  we  deduce  the  following 
eigenvalue  problem  expressed  in  terms  of  the  electric  field 

rot  (  rot  e)-€^€QM.QW^e  =  0  infi  (5) 

n  /V  e  =  0  on  r  (6) 

A  weak  formulation  of  (5) -(6)  holds  in  Eg  [11]: 

Jq  rot  e.rot  e'dfi  -k^Jj^c^e.e'dfi  =  0  \/e'GEg(7) 

where  ^ 

Eg  =  ^  e  e  [I?  rot  e  G  [L^  (fi)  ]  ,  n  -  e  Ir  =  0 

and  where  the  values  of  k  (k^=(o^€gM-o) 
wavenumbers . 

SPURIOUS  MODES 

The  searched  resonant  fields  (u  >  0  )  are 

theoretically  divergence  free  since  from  (5)  it  follows: 

div  €  ^e  =  -  div  (  rot  (rot  e) )  =  0  (8) 

}c 

Moreover,  for  a  simply-  connected  region  fi  ,  the  only  field 
corresponding  to  o)=0  (static  field)  satisfying  eguations 
(1) , (3)  and  the  condition  div  €^e=0  is  e-0. 

The  trouble  arises  when  discretizing  (7)  with 
classical  finite  elements  (for  example  nodal  vector 
elements)  :  a  matrix  with  many  eigenvalues  being  zero  is 
obtained  (0  is  a  highly  degenerate  eigenvalue)  [3]  [5].  The 
numerical  approximations  of  this  value  k^=0  are  difficult 
to  isolate  from  the  meaningful  lowest  non-zero  eigenvalues 
(  k^  >  0) ;  especially  when  the  number  of  degrees  of  freedom 
increases.  Most  of  them  do  not  satisfy  div  €^e=0  and  then 
are  unacceptable  as  solutions  of  Maxwell's  equations.  The 
resulting  set  of  solutions  is  a  mixture  of  physical  modes 
and  numerical  spurious  ones. 

With  edge  elements,  as  we  shall  see  in  the  next 
section  such  a  situation  doesn't  occur  . 

FINITE  ELEMENT  DISCRETIZATION 

Mixed  finite  elements  [7]  are  used  for  the  numerical 
approximation  of  (7) .  Let  (  with  boundary  )  be  the 

discretization  of  with  tetrahedra.  The  "edge  elements” 
have  the  following  properties: 
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The  degrees  of  freedom  e^  and  the  trial  functions 
are  associated  with  the  mesh  edges.  For  every  edge  "a" 
containing  the  nodes  “i''  and  "j"  : 

.  e^  is  the  circulation  of  e  along  "a”. 

®a  =  Ja 

.  Wg  can  be  expressed  in  terms  of  the  barycentric 
functions  X,-  and  Aj  as  : 

Wg  =  Aj  gradAj  -  Aj  gradA,-  (10) 

In  each  tetrahedron  e  is  : 

e(r)  =  a  +  (3  ^  r  (11) 

where  a  and  p  are.  three-components  constant  vectors  and  r 
is  the  vector  (x,y,z). 


We  introduce  the  space  : 


E 


Oh 


(12) 


For  every  e  in  Eq^  the  tangential  part  of  e  is  continuous 
across  tetrahedra  interfaces.  The  approached  problem  is  to 
find  e  in  Eg^  so  that  : 

rot  e.rot  e'dfi^  -  m  ^  ^e.e'dfi^  =0  V  e'  G  Eg^  (13) 

h  V 

Finally,  we  have  to  solve  a  generalized  algebraic 
eigenvalue  problem  of  the  form  : 

A  u  =  B  u  (14) 

A  ("stiffness  matrix")  and  B  ("mass  matrix")  have 
dimensions  n^x  n^  where  n^  is  the  number  of  edges  in  the 
finite  element  mesh. 


Remark:  with  "edge  elements"  all  the  numerical  solutions 
corresponding  to  k^  >0  are  "weakly  divergence  free"  :  no 
spurious  mode  occur.  The  reason  is  the  following  : 

Let  <p  be  any  linear  combination  of  the  barycentric 
functions  A,.  (  <p  piecewise  affine  )  and  <p  =0  on  F;  it  can 
be  showed  that  all  the  fields  b'  of  the  form  e'=  grad  <p 
are  in  Egj,  [11]  [12].  Then  they  can  be  chosen  as  admissible 
test  fields  in  (13).  Rewriting  e'=  grad<p  in  (13)  leads  to  : 
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€pe.  grad<p  =  0  (15) 

Equation  (15)  is  equivalent  to  say  that  for  every  k^>  0  the 
solution  div  €^e  =  0  is  verified  in  the  distribution  sense. 
An  integration  by  parts  of  (15)  shows  that  for  every  inner 
node  "i”  of  ,  an  integral  involving  the  sum  of  jumps  of 
e.n  through  facets  of  tetrahedra  having  node  '‘i”  in 
common  is  zero.  This  property  is  important  since  it  gives 
an  average  local  divergence  condition.  Such  result  does  not 
exist  in  case  of  classical  finite  elements  (nodal  ones) 
since  the  fields  of  the  form  e'  =  grad<p  are  not  in  the 
space  of  the  test  functions  E^^  . 

NUMERICAL  RESULTS 


a.  Empty  rectancfular  cavity. 

An  empty  rectangular  cavity  with  perfectly  conducting 
walls  was  modelled  with  the  above  developed  technique.  The 
cavity  has  dimensions:  a=  0.4m,b  =0.3m,c=l.m 
(figure  1-a) .  A  quarter  of  the  cavity  (figure  1-b)  was 
analyzed  using  220  tetrahedra.  Symmetry  conditions  were 
prescribed  on  faces  x  =  0  and  z  =  0  . 


a) 


b) 


Figure  1  Studied  air-filled  cavity 


The  problem  is  symmetrical  in  x,  y,  z;  so  the  fields 
can  be  expressed  as  TE  (transverse  electric)  or  TM 
(transverse  magnetic)  to  any  one  of  these  coordinates  [13]. 
It  is  conventional  to  choose  the  longer  dimension  along  the 
z  direction.  Analytical  solutions  are  then  labelled  TE^ 
(modes  whose  electric  field  has  no  z-component  )  and  TM^ 
(modes  whose  magnetic  field  has  no  z-component) .  The 
corresponding  resonant  wave  numbers  are  : 


137 


,  5 

)c2  =  it2  (  —  +  _  +  ±L  ) 
a2  b2  c2 

The  algorithm  for  solving  the  matricial  eigenvalue 
problem  is  based  on  the  classical  QR  method. 

Computations  were  performed  on  a  DN4000  Apollo 
workstation;  about  10  minutes  are  necessary  to  solve  the 
entire  problem.  The  six  lowest  computed  modes  are  shown  on 
Table  1. 


Mode 

Wavenumber 
k  computed 

Wavenumber 
k  analytical 

Error 

(%) 

TEioi 

8.365 

8.458 

1. 

TE,„  (TM,„  ) 

13.243 

13.461 

1.6 

13.488 

13.461 

0.2 

™,05 

12.327 

12.268 

TE,„  (TM,,,  ) 

16.273 

16.129 

0.9 

16.356 

16.129 

1.4 

Table  1 

Nimerical  Results 


This  simple  case  is  known  to  give  spurious  solutions 
when  solved  with  classical  finite  elements  [2]  [14].  Here 
no  spurious  mode  is  observed.  Moreover  the  relative  error 
never  exceeds  2%. 


Figure  2  Figure  3 

Electric  field  (TE^q,  mode  Electric  field  (TE,,,  mode 

in  (x,y)  plane  for  z=0.5  m)  in  (x,y)  plane  for  z=0.5  m) 

Figures  2  and  3  show,  in  the  quarter  of  the  cavity, 
the  distribution  of  vector  fields  for  the  TE^q^  mode  and 
the  TE^,,  mode  respectively.  The  plane  of  symmetry  is  on 
the  right  of  the  figures.  On  this  plane  e  is  tangential  and 
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on  the  others  e  is  normal.  In  each  tetrahedron  e  is 
represented  with  an  arrow  in  the  centre  of  gravity  and  the 
length  of  the  arrow  is  proportional  to  the  magnitude. 

b. Inhomogeneous  dielectric  loaded  cavities 

Two  examples  of  inhomogeneous  loaded  cavities  were 
analyzed. 

The  first  one  (  figure  4  )  is  the  preceding  cavity 
with  dielectric  discontinuities  in  one  direction  only.  The 
relative  permittivity  of  the  dielectric  material  is 

=16  .  For  a  quarter  of  the  cavity  230  tetrahedra  were 
used.  The  theoretical  eigenvalue  for  the  dominant  mode 
(lowest  eigenvalue)  is  known  [15]  :  k  =  2.5829.  An  error  of 
0.4%  was  found. 


‘  Figure  4 

Cavity  with  dielectric  block  (example  1) 

The  second  cavity  (figure  5-a)  is  the  cavity  of  a. 
with  dielectric  discontinuities  in  three  dimensions  (  = 
16) .  The  quarter  (figure  5-b)  was  modelled  with  240 
tetrahedra  (figure  5-c) .  Computing  time  is  about  15 
minutes . 


0.25  m  a)  b) 
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C) 

Figure  5 

Cavity  with  dielectric  block  (example  2) 


No  analytical  result  is  available  but  comparison  with 
already  published  computed  values  is  possible;  results  for 
the  dominant  mode  are  the  following; 

Source  k  (computed  ) 


Ref  [1]  5.60 

Ref  [15]  5.529 

Ref  [16]  4.907 

Presented  Method  5.102 


All  these  results  agree  within  roughly  10%.  Some 
others  structures  should  be  modelled  in  order  to  make  a 
comparison  more  satisfactory  between  all  the  methods. 
However  the  mixed  finite  elements  used  here  are  well 
adapted  in  case  of  dielectric  materials  because  they  imply 
the  tangential  continuity  of  the  electric  field  across 
interfaces  and  take  account  of  the  discontinuity  of  the 
normal  component  . 


Figure  6 

Electric  field  (dominant  mode) 

In  the  (x,y)  plane  In  the  (y^z)  plane 

for  2=0.5  m  for  x=0.2  m 
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No  spurious  mode  occured;  all  the  computed  field  correspond 
to  resonant  fields.  Figure  6  shows  electric  vector  fields  in  the 
quarter  of  the  cavity. 

CONCLUSION 

"Edge  elements”  (a  class  of  mixed  finite  element)  have  been 
used  to  model  empty  and  dielectric  loaded  cavities.  The  first 
resonant  frequencies  were  computed;  comparison  with  analytical 
values  or  results  published  in  previous  papers  shows  the 
efficiency  of  the  method. 

These  elements  avoid  all  the  well  known  "spurious  modes"  and 
seem  very  promising  for  the  study  of  more  complicated  problems  in 
high  frequency  applications 
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ABSTRACT 

The  fast  concurrent  implementation  of  a  FORTRAN  method  of  moments  analysis  of  the 
electromagnetic  properties  of  an  array  of  tapered  slot  antennae  is  discussed.  Decomposition  of  an 
existing  FORTRAN  algorithm  for  calculation  of  the  currents  induced  by  an  incident  radiation  field  in 
an  infinite  array  of  tapered  slot  antennae  is  described.  The  problem  was  distributed  across  an  array  of 
INMOS  transputers,  yielding  significant  speed-up  over  a  single  CPU.  This  decomposition  was 
relatively  simple  to  implement,  can  readily  be  scaled  to  larger  processor  arrays  virtually  indefinitely , 
and  promises  linear  speed-up  with  the  number  of  processors  in  the  array. 

1 .  Introduction 

A  growing  number  of  electromagnetic  analysis  problems  are  being  formulated  for  processing  on 
parallel  processing  computers,  and  transputer  arrays  represent  one  of  the  least  costly  parallel  computer 
systems  available  to  the  EM  analyst.  Recently,  tutorial  papers  have  appeared  to  describe  the  successes 
and  difficulties  one  may  encounter  in  solving  various  types  of  problems  [1,  2].  The  objective  of  this 
paper  is  to  describe  the  way  in  which  a  FORTRAN  moment  method  analysis  was  converted  from  a 
code  for  a  single  CPU  to  a  code  that  utilises  several  cpus  with  very  high  efficiency. 

The  particular  analysis  that  was  converted  is  for  infmite  arrays  of  endfire,  tapered  slot  antennas  [3]-[5]. 
The  currents  flowing  on  the  metallic  fins  that  comprise  an  infinite  array  of  tapered  slot  antennae  (figure 
1)  are  determined  by  using  the  method  of  moments  to  solve  the  electric  field  integral  equation. 
Floquet's  theorem  is  used  so  that  only  a  unit  cell  of  the  structure  must  be  considered,  and  the  version 
of  the  analysis  being  performed  here  uses  piecewise  sinusoidal  rooftop  basis  and  testing  functions. 
Most  (over  95%)  of  the  time  taken  to  calculate  the  currents  induced  in  the  antennae  is  spent  filling  the 
impedance  matrix  for  the  50  to  200  unknowns  that  are  used  to  model  the  current.  Since  this  particular 
analysis  is  not  the  main  topic  of  the  paper,  only  a  few  sample  results  are  included. 


y4- 
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Figure  1  Tapered  slot  antennae 

Use  of  the  numerical  model  for  antenna  design  and  development  requires  computing  the  antenna 
currents  for  many  scan  angles  and  frequencies  of  operation.  Computations  for  each  set  of  parameters 
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are  lengthy  but  tractable  with  modest  CPUs;  typical  cases  require  several  minutes  to  a  few  hours  on  a 
25-MHz  80386/80387  PC.  However,  typical  applications  require  that  dozens,  or  perhaps  hundreds, 
of  parameter  sets  be  analyzed.  It  is  the  objective  of  this  paper  of  demonstrate  that  EM  analysts  can 
readily  decompose  a  FORTRAN  algorithm  of  this  type  to  be  distributed  to  an  economical  ttansputer 
array.  The  decomposition  maintains  the  input  and  output  features  familiar  to  the  user  while  distributing 
the  computations  for  each  parameter  set  to  a  different  processor  via  a  "general  purpose" 
communications  harness  that  is  easily  adapted  to  different  FORTRAN  algorithms.  A  request-driven 
transputer  farming  paradigm  is  used,  which  provides  efficient  utilization  of  the  transputer  array  for  this 
type  of  computation.  Thus,  we  demonstrate  how  easily  the  total  processing  speed  of  the  transputer 
array  can  be  used  to  the  benefit  of  EM  analysis  codes  that  were  not  originally  intended  for  parallel 
processing  as  well  as  to  the  benefit  of  new  codes  that  may  be  written  specifically  to  exploit  the 
parallelism  in  computations  for  a  single  parameter  set.  The  conttoller/worker  model  is  described  in 
detail  in  the  implementation  section  below. 

2 .  Parallel  Implementation 

2.1.  Paradigms  for  Parallelism 

Parallel  systems  may  be  described  at  various  levels  of  abstraction,  and  it  proves  convenient  to  use  the 
subclasses  conceptual,  physical,  and  logical  architectures  to  describe  more  accurately  what  is  meant  by 
'parallel  architectures'.  This  is  more  than  mere  taxonomy;  understanding  parallel  processing  and 
effectively  managing  parallel  systems  design  and  implementation  requires  different  views  to  be  taken 
for  different  purposes,  and  thinking  in  terms  of  these  architectural  subclasses  has  been  found  very 
useful  in  practice. 

Conceptual  architectures  involve  the  most  abstract  aspects  of  parallel  processing  systems,  in  particular 
the  type  of  computation  being  undertaken.  Most  conceptual  architectures  fall  into  one  of  four  classes: 
control-driven,  data-flow,  object-based,  and  logic-based. 

Physical  architectures  deal  with  the  actual  physical  implementation  of  such  systems.  The  categorisation 
of  physical  architectures  preferred  is  that  introduced  by  Hynn  [6].  It  is  a  classification  which 
distinguishes  four  basic  combinations  of  instruction  streams  and  data  streams  (SISD,  SIMD,  MISD, 
and  MIMD).  Flynn's  scheme  is  simpler  than  most,  and  is  very  useful  for  classifying  machines  which 
are  clearly  different.  This  work  uses  a  machine  in  the  last  of  the  above  four  categories  -  a  transputer 
array. 

The  logical  architecture  of  a  system  is  the  architecture  seen  by  applications  software  designers  and 
programmers.  There  are  three  main  types  of  logical  architecture  appropriate  to  MIMD  machines  - 
farming,  geometric  and  algorithmic  [7,8]. 

Farming  is  a  controller/worker  paradigm  in  which  an  overall  problem  is  split  into  conveniently  sized 
work  packets  which  are  disuibuted  over  a  number  of  worker  processors  in  such  a  way  as  to 
dynamically  balance  the  work  load  across  the  processors  in  the  farm.  The  controller  knows  about  the 
work  that  has  to  be  done,  and  the  worker  knows  how  to  do  it  -  each  worker  has  an  identical  copy  of 
the  code.  The  controller  sends  batches  of  work  to  each  worker,  the  latter  executes  its  algorithm  on  a 
batch  of  data  then  returns  the  results,  starting  to  work  on  the  next  batch  pending  any  message 
transmission.  The  primary  advantage  of  the  task  farm  is  the  fact  that  it  dynamically  balances  the 
loading  of  an  overall  problem  across  a  network  of  processors.  A  potential  disadvantage  is  that  each 
worker  has  to  have  a  copy  of  the  complete  code,  which  can  lead  to  problems  if  memory  is  limited. 

An  architecture  of  this  type  is  most  appropriate  to  an  application  in  which  firstly,  the  work  packets  may 
require  different  amounts  of  processing  -  that  is  to  say  the  amount  of  work  for  each  is  packet- 
dependent.  Secondly,  and  perhaps  more  obviously,  the  overall  task  should  naturally  lend  itself  to 
temporal  concurrency.  This  means  that  as  far  as  is  possible,  the  processing  of  any  individual  packet 
should  not  depend  upon  the  processing  of  any  other  packet.  Complete  independence  may  not  actually 
be  possible;  packets  may  have  some  positional  or  other  relationship  to  each  other.  However,  such  an 
architecture  does  demand  that  the  final  result  does  not  depend  upon  the  order  in  which  the  packets  are 
processed. 
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In  contrast  with  the  automatic  run-time  load-balancing  in  a  farm,  geometric  and  algorithmic 
architectures,  which  exploit  data  concurrency  and  algorithm  concurrency  respectively,  have  to  be  load- 
balanced  at  design  time. 

In  geometric  parallelism,  the  parallelism  inherent  in  the  data  is  exploited.  The  problem  will  generally 
have  an  underlying  geometrical  structure  (this  is  often  very  similar  to  SIMD  array  processing),  and 
when  this  is  the  case,  the  array  topology  reflects  the  data  topology.  In  geometric  parallelism,  as  in 
farming,  each  processor  has  an  identical  copy  of  the  complete  code  for  the  whole  program,  but  each 
works  on  separate  pre-defined  portions  of  the  data  space,  i.e.,  each  processor  is  responsible  for  a 
specific  area  or  volume  of  the  data  set. 

Algorithmic  parallelism  is  the  case  in  which  the  parallelism  inherent  in  the  algorithm  is  exploited,  e.g. 
with  a  pipeline.  In  a  distributed  memory  machine,  each  processor  has  its  own  segment  of  the  code,  and 
this  is  different  from  processor  to  processor.  The  topology  of  the  machine  typically  reflects  the 
topology  of  the  data  flow  graph  on  which  the  algorithm  is  based. 

Many  architectures  are,  in  fact,  hybrids  of  two  or  all  three  types,  but  our  decomposition  is  purely  a 
farm.  The  logical  architecture  is  described  in  more  detail  a  little  later.  Before  giving  the  description  we 
present  a  brief  summaiy  of  the  various  techniques  used  in  implementation  of  parallel  software. 

2.2.  Implementing  Parallel  Software 

There  are  several  distinct  approaches  to  implementing  parallel  software.  The  quickest,  and  the  one  used 
in  this  case,  is  to  use  existing  codes  within  a  specially  created  communication  harness.  Possible 
parallelism  is  identified  at  the  procedure  level,  separate  code  segments  are  isolated,  and  then  a 
communication  harness  is  constructed.  Slight  modifications  may  be  needed  to  the  existing  code  -  the 
extent  of  this  depends  on  the  manner  in  which  the  original  code  was  written.  According  to  how  well 
the  existing  code  is  written,  this  approach  usually  requires  less  effort  than  either  implementing  the 
same  algorithm  in  a  more  appropriate  language  (e.g.,  Occam)  or  than  making  a  completely  fresh  start 
with  respect  to  algorithm  and  language,  and  provides  a  significant  return. 

2.3.  Transputers  and  Parallelism 

Confusion  is  often  caused  in  discussions  of  transputer-based  software  by  the  ambiguous  use  of  the 
term  parallelism.  Three  distinct  forms  of  parallelism  may  be  observed  in  a  multi-transputer 
implementation. 

Firstly  there  is  the  case  of  separate  transputers  running  separate  processes.  This  is  genuine 
simultaneous  operation  of  concurrent  or  parallel  processes,  similar  to  having  two  or  more  independent 
computers  running  side  by  side. 

Secondly,  a  single  transputer  may  be  running  several  processes  in  a  time-sliced  mode.  This  is 
analogous  to  the  forms  of  pseudo-parallelism  performed  on  many  conventional  machines,  e.g.,  those 
running  processes  under  an  operating  system  such  as  Unix  or  VAX  VMS.  Time  slicing  is  performed 
on  a  time  scale  which  is  short  compared  to  the  process  lifetimes.  It  should  be  remembered  that  each  of 
the  parallel  processes  on  a  single  transputer  is  competing  for  epu  time,  and  in  the  same  way  that 
another  user  on  a  VAX  makes  operation  seem  slower,  every  extra  process  will  slow  down  the 
operation  of  those  already  there. 

The  time-slice  mechanism  of  the  transputer  is  implemented  in  hardware,  and  is  particularly  efficient 
compared  to  many  similar  systems.  In  particular,  processes  waiting  for  communication  are  removed 
from  the  active  process  list  (descheduled)  and  incur  no  epu  overheads. 

The  third  form  of  parallelism  concerns  the  hardware  of  the  transputer.  The  TSOO  contains  a  fixed-point 
epu,  a  floating-point  processor,  and  four  bidirectional  serial  communications  link  engines;  all  of  these 
elements  can  function  simultaneously  and  independently.  This  is  particularly  important  in  relation  to 
communications.  If  the  time  required  to  pass  a  message  is  less  than  the  processing  time  required  to  do 
the  work  generated  by  the  previous  message,  then  the  communications  can  be  done  in  the  background 
with  little  impact  on  the  time  required  to  generate  the  desired  results. 
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The  link  engines  operate  by  DMA.  For  a  message  to  pass  between  two  processors,  the  cpu  must 
provide  the  link  engine  with  a  start  address  and  a  byte  count.  This  cpu  overhead  for  communications 
means  that  it  is  more  efficient  to  pass  long  messages  (e.g.,  a  few  kbytes)  than  short  ones  (e.g.,  single 
bytes). 

In  this  paper  we  are  primarily  concerned  with  the  exploitation  of  the  genuine  parallelism  of  running 
separate  processes  on  separate  processors. 

3.  Algorithm  FOR  THE  CPU 

The  analysis  was  originally  coded  in  FORTRAN  for  implementation  on  a  single  CPU.  The  code  is 
comprised  of  a  main  program  which  performs  data  input  and  output,  and  controls  the  computations  by 
calling  subroutines  that  fill  the  impedance  matrix  and  solve  the  system  of  equations.  A  typical  matrix 
element  is  obtained  by  evaluating  an  expression  of  the  form  [9]: 


(nAVm)  sin(pmna) 
(Pmna)  [cos(Uoa)  -  cos  pmna] 


sin(VniWq)  cos(nAZq)[cos(nAhq)  -  cos(khq)] 

Vm(k2-n2A2)sin(k  hq) 


eJVmyq 


^  sin(nAzD)sin(nAWp)[cos(VinhD)  -  cos  (k  hp)]  ^y^yp 


(nA)(k2-VM  sin(k  hp) 


k  =  (0 

00.  =  main  beam  direction 

a,b  =  array  grid  spacings  in  x  and  y  directions 
A  =  increment  of  the  spectral  variable 
yq,  Zq  =  coordinates  of  center  of  basis  function  q 
yp,Zp  =  coordinates  of  the  center  of  testing  function  p 
hq,Wq  =  half-length  and  half- width  of  basis  function  q 
hp,Wp  =  half-length  and  half-width  of  testing  function  p 

The  upper  limits  of  the  summations  required  for  convergence  are  typically  +/-  15  for  m  and  500-1500 
for  n.  The  large  range  of  the  upper  limit  for  n  is  due  to  this  parameter  being  sensitive  to  the  antenna 
geomet^.  The  FORTRAN  code  uses  the  common  practice  of  storing  various  terms  that  comprise 
Kmn.  3nd  T^^  so  that  they  are  not  recomputed  when  needed  for  successive  sets  of  the  indices. 
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Total  storage  requirements  are  kept  to  a  reasonable  level  (1-3  MBytes)  by  proper  nesting  of  the  loops 
on  m,  n,  p  and  q.  The  problem  possesses  an  inherent  parallelism  because  the  computations  are  usually 
perfomied  for  several  scan  angles,  which  provides  a  simple  way  to  divide  it  up  into  independent  sub¬ 
problems.  Except  for  a  few  initialising  calculations,  the  various  scan  angles  do  not  share  any  common 
data,  but  the  code  that  is  executed  for  each  angle  is  identical. 

4,  The  fortran  farm 

4.1.  Toplevel  architecture 

The  architecture  used  in  this  case  is  shown  in  figure  2  -  a  ring.  Data  and  work  packets  pass  round  the 
farm  in  the  same  direction.  All  file  access  and  control  is  provided  by  the  driver  and  the  calculations  are 
performed  by  the  worker  processors.  The  communications  harness  was  imtten  in  the  parallel 
processing  language  Occam,  and  all  the  original  file  access,  control  and  calculation  code  was  retained. 
The  result  is  about  600-700  lines  of  Occam  and  about  2000  lines  of  FORTRAN.  Occam  was  chosen 
for  the  communications  harness  because  it  is  a  natural  language  for  describing  parallel  communications 
within  a  computer  program  and  it  was  the  quickest  means  available  at  the  time  of  creating  an  efficient 
communications  harness.  Further  data  on  the  transputer  system  are  contained  in  table  1. 


Table  1  System  characteristics 

Hardware  Cotnpligis 

Meiko  MIO  computing  surface  Inmos  Occam  compiler  for  the 

communications  harness 

T800  transputers  (20  MHz,  4  cycle  RAM) 

T800  maximum  computing  power  =  IMflop  u 

VAX3800  maximum  computing  power  *=  3-4  Mflops  Meiko  Fortran  77  compiler  for  the 

main  program 


results 


Fi^re  2  Toplevel  architecture 

In  this  case  the  most  convenient  decomposition  was  to  perform  the  calculations  for  different  scan 
angles  on  different  processors.  This  is  the  simplest  approach  possible,  and  is  also  the  one  which  yields 
the  highest  ratio  of  computation  load  to  communications  load.  The  work  packets  passed  out  to  the  farm 
consist  of  little  more  than  the  scan  angle  at  which  the  calculations  are  to  be  performed  -  a  few  bytes  of 
data  at  most,  whereas  the  calculations  for  each  angle  require  between  one  and  two  hours  of  T800  cpu 
time. 


4.2.  Communications  paradigm 

The  simplest  type  of  communications  harness  for  a  farming  decomposition  is  one  in  which  the  driver 
processor  passes  work  packets  out  to  the  farm  as  quickly  as  possible,  and  each  worker  processor 
accepts  data  if  it  is  free  to  do  some  work,  or  passes  it  on  if  not.  This  flood-fill  approach  to  supplying 
the  workers  with  data  for  processing  is  not  appropriate  in  this  case,  because  it  is  necessary  for  each 
worker  to  have  buffer  processes  in  which  work  packets  must  be  stored  temporarily  whilst  the 
processor  decides  whether  to  pass  on  the  data  or  accept  it  for  processing  -  the  minimum  possible 
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number  of  buffers  is  one.  This  is  not  usually  a  problem  for  decompositions  in  which  the  number  of 
processors  in  the  farm  is  one  or  two  orders  of  magnitude  less  than  the  total  number  of  work  packets. 
However,  in  a  case  where  the  total  number  of  work  packets  to  be  processed  is  only  twice  or  three 
times  the  number  of  processors  in  the  farm,  the  flood-fill  approach  can  have  very  serious 
disadvantages.  Work  packets  are  often  left  stranded  in  buffers  in  such  a  way  that  at  the  end  of  a 
processing  run  the  last  two  or  three  packets  are  processed  one  after  the  other  by  the  same  processor. 
This  is  not  noticeable  in  cases  where  there  are  hundreds,  or  even  thousands,  of  individual  work 
packets,  each  of  which  takes  a  few  milliseconds  to  process.  If  there  are  only  a  few  tens  of  packets  and 
a  similar  number  of  processors,  as  there  are  in  this  particular  decomposition,  then  the  tot^  compute 
time  can  be  significantly  greater  than  if  the  farm  were  used  with  the  maximum  possible  efficiency. 

We  have  adopted  an  approach  which  removes  the  necessity  for  any  buffering  on  the  worker  processors 
-  the  farm  is  request  driven.  The  driver  only  passes  work  packets  out  to  the  workers  when  they 
specifically  request  it.  Each  packet  is  tagged  with  the  address  of  the  requesting  worker  and  no  other 
processor  may  accept  that  packet.  This  approach  can  reduce  farm  efficiency  if  the  work  packets  are 
large  because  there  is  then  a  significant  delay  between  a  worker  sending  out  a  request  and  receiving  its 
work  package,  during  which  time  the  processor  stands  idle;  but  for  this  problem  the  work  packets  are 
at  most  ten  bytes,  and  for  a  transmission  time  around  ten  microseconds  per  processor  the 
communications  delay  is  insignificant. 

4.3.  Performance 

Typical  processing  runs  were  over  20  to  55  scan  angles  with  10  to  20  worker  processors  available. 
The  results  below  are  for  a  particular  example  of  52  scan  angles  (0  =  0®  to  85®  in  5®  steps  and  <|)  —  0° 
to  90°  in  45®  steps)  distributed  over  17  processors:  the  time  taken  for  the  processing  run  was  7.61 
hours.  The  same  processing  run  on  a  VAX  3800  with  a  floating  point  accelerator  took  37.32  hours, 
even  without  the  benefit  of  the  extremely  efficient  FORTRAN  compiler  available  on  the  VAX,  the 
transputer  system  is  four  times  as  fast.  This  increase  in  compute  power  allowed  overnight  production 
of  results  which  would  otherwise  have  taken  two  days  of  dedicated  VAX  3800  cpu  time. 

The  VAX  3800  is  chosen  for  comparison  because  its  characteristics  are  widely  publicized  and  these 
machines,  or  similar  ones,  are  often  available  to  EM  analysts.  Obviously,  newer  workstations  could 
outperform  the  VAX  3800,  but  newer  transputers  also  could  ouqierform  the  T800. 


4.4.  Sample  Computed  Results 

Typical  results  for  the  input  impedance  of  the  antenna  depicted  in  figure  3  are  shown  in  figures  4  and 
5.  These  curves  are  similar  to  those  one  might  obtain  for  any  well-behaved  phased  array  antenna.  The 
resistance  versus  scan  is  reasonably  constant  near  broadside  and  then  drops  to  zero  at  endfire.  In  the 
H-plane  ((j)  =  0)  and  the  E-plane  (<1)  =  90),  the  reactance  changes  are  opposite.  The  same  array, 
operating  at  3.725  GHz,  exhibits  at  grating  lobe  at  54.9®  in  the  principal  planes.  The  reflection 
coefficient  for  this  case  is  shown  in  Figure  6.  The  data  in  Figures  4,  5  and  6  were  computed  by  using 
77  rooftop-like  current  modes  on  the  metal  fin.  These  modes  and  a  closely  related  formulation  of  the 
electromagnetic  problem  are  described  in  [9],  as  are  issues  related  to  convergence  and  validation  of  the 
analysis.  To  summarize  those  findings,  the  summation  limits  mentioned  above  are  adequate  for 
computing  the  impedance  to  an  accuracy  that  agrees  within  1  or  2  ohms  with  published  data  for  test 
cases,  such  as  dipoles  and  monopoles,  and  also  agrees  well  with  waveguide  simulators  for  other 

cases. 

5.  CONCLUSIONS 

We  calculate  that  the  ratio  of  communications  to  compute  load  for  this  particular  application  is  so  small 
(a  few  milliseconds  versus  half  an  hour)  that  the  farm  size  could  easily  be  increased  to  several  hundr^ 
processors  without  significant  (or  even  visible)  degradation  of  efficiency.  This  implies  speed¬ 
up.  This  assumes  either  that  a  more  finely  grained  decomposition  for  this  algorithm  could  w  found 
without  excessively  increasing  the  communications  load  on  the  farm,  or  that  one  would  desire  to 
compute  the  antenna  response  either  for  many  thousands  of  scan  angles,  of  for  several  tens  of  scan 
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angles  at  many  different  frequencies.  Increasing  the  number  of  processors  is  easy,  finances  permitting 
-  in  the  software,  a  single  constant  defining  the  number  of  workers  is  changed  and  the  communications 
harness  is  recompiled  -  the  FORTRAN  is  completely  unaffected.  We  estimate  that  a  32  processor 
system,  offering  an  order  of  magnitude  speed-up  over  a  VAX  3800,  complete  with  all  necessary 
software  and  additional  hardware  could  be  purchased  for  less  than  half  the  price  of  the  VAX. 

An  alternative  to  increasing  the  number  of  processors,  or  redesigning  and  recoding  the  algorithm  to 
suit  a  finer  grained  decomposition,  is  to  use  different  hardware.  Examples  of  high  performance  cpus 
available,  or  soon  to  be  available,  are  the  INTEL  i860  and  the  INMOS  T9000  transputer.  These  would 
both  offer  performance  improvements  over  the  T800  of  at  least  a  factor  of  ten  -  the  i860  is  already 
available  (with  communications  handled  by  transputers),  and  the  T9000  is  expiected  to  be  available  late 
in  1991.  The  i860-transputer  hardware  can  be  obtained  for  around  £10,000  per  i860:  as  yet  the  truly 
exceptional  performance  of  this  chip  (upto  60Mflops  -  60  x  T800  transputer  power)  is  only  available  to 
those  who  are  willing  to  hand  craft  their  algorithms  in  the  appropriate  assembler.  Currently  available 
compilers  achieve  between  5  and  10  Mflops.  Between  the  times  of  writing  the  original  manuscript  and 
the  revisions,  the  cost  of  i860  boards  has  been  reduced  more  than  50  percent.  However,  at  present,  it 
appears  that  the  30  Mflops/£20,000  cost  effectiveness  of  the  transputer  farm  employing  TSOO's  is 
superior  to  that  of  multiple  i860's  that  would  be  required  to  achieve  30  Mflops  with  currently  available 
compilers.  This  means  that  the  i860  is,  as  yet,  no  more  cost  effective  than  the  transputer  -  we  are  not 
yet  able  to  comment  on  the  T90(X). 
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